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Overview 


This  project  has  received  support  from  ARO  to  advance  the  state  of  the  art  in  the  CFD 
solution  of  the  Euler  and  Navier  Stokes  equations,  with  the  specific  objective  of 
increasing  accuracy  in  the  computational  prediction  of  compressible  flows  of  Army 
interest.  To  reach  this  objective,  this  project  has  developed  an  innovative  implicit 
finite  element  CFD  algorithm  that  relies  on  multi-dimensional  characteristics  theory. 
By  firmly  resting  on  the  physics  and  mathematics  and  acoustics  and  convection,  the 
wave  propagation  mechanisms  within  fluid  flows,  the  algorithm  is  designed  to 
generate  accurate  solutions  and  minimize  induced  artificial  diffusion. 

During  this  project,  a  40-page  paper  detailing  project  results  and  entitled:  "A  CFD  Euler 
Solver  from  a  Physical  Acoustics-Convection  Flux  Jacobian  Decomposition"  has  been 
published  in  the  International  Journal  for  Numerical  Methods  in  Fluids.  This  paper  is 
included  as  Appendix  G  of  this  report.  Furthermore,  the  principal  investigator  has 
signed  a  publishing  contract  with  John  Wiley  &  Sons  to  write  an  advanced  CFD  book 
that  will  feature  the  results  obtained  in  this  project. 
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SUMMARY 

This  report  details  for  the  3-D  Euler  and  Navier-Stokes  equations  an  intrinsically  infinite  directional  up¬ 
stream  formulation  that  rests  on  the  mathematics  and  physics  of  multi-dimensional  acoustics  and  convec¬ 
tion.  Based  upon  characteristic  velocities,  this  formulation  introduces  the  upstream  bias  directly  at  the 
differential  equation  level,  before  the  spatial  discretization,  within  a  characteristics-bias  governing  system 
A  conventional  centered  discretization  of  this  system  on  given  grids  directly  yields  an  optimal  discretely 
conservative  and  multi-dimensional  upstream  approximation  for  the  Euler  equations.  Through  a  decom¬ 
position  of  the  Euler  flux  divergence  into  multidimensional  acoustics  and  convection  acoustics  components, 
the  formulation  induces  consistent  upstream  bias  along  all  directions  of  spatial  wave  propagation,  with 
anisotropic  variable-strength  upstreaming  that  correlates  with  the  spatial  distribution  of  characteristic  ve¬ 
locities.  With  the  objective  of  minimizing  induced  artificial  diffusion,  the  formulation  non-linearly  induces 
upstream-bias  essentially  locally  in  regions  of  solution  discontinuities,  whereas  it  decreases  the  upstream-bias 
in  regions  of  solution  smoothness.  The  discrete  equations  originate  from  a  finite  element  discretization  of  the 
characteristic-bias  system  and  are  integrated  in  time  within  a  compact  block  tri-diagonal  matrix  statement 
by  way  of  an  implicit  non-linearly  stable  Runge-Kutta  algorithm  for  stiff  systems. 

KEY  WORDS:  CFD;  characteristics;  multi-dimensional  upwind;  implicit;  finite  elements; 


1  Introduction 

The  CFD  community  is  still  eagerly  pursuing  a  multi-dimensional  Euler  solver  for  computing  real¬ 
istic  flows  on  arbitrary  grids.  Many  flnite  element,  difference  and  volume  algorithms  have  remained 
somewhat  independent  from  the  physics  of  acoustics  and  convection,  the  wave  propagation  mech¬ 
anisms  within  gas  dynamic  flows.  The  dissipation  mechanisms,  or  upwind  schemes,  within  these 
algorithms  have  been  developed  at  the  discrete  level,  in  connection  with  a  specific  grid  or  pattern 
of  computational  cells. 

Several  finite  element  solvers  have  either  utilized  modifications  of  the  test  space  or  introduced 
Taylor’s  series  based  dissipation  terms,  e.g.  TWS,^  to  generate  stable  algorithms.  The  mathemat¬ 
ical  developments  in  these  fundamental  contributions  have  remained  independent  of  characteristics 
theory.  Upwind  finite  element  methods  for  scalar  equations  have  also  been  developed,^  including 
the  Streamline  Upwind  Petrov-Galerkin  (  SUPG  )  formulation^"^  for  also  known  as  the  Streamline 
Diffusion  (  SD  )  method.  Extensions  to  systems  are  recognized  to  remain  heuristic,®  the  induced 
upwinding  is  not  necessarily  in  the  streamline  direction,  and  additional  “shock  capturing”  terms 
are  needed  for  computing  essentially  non-oscillatory  shocked  flows. 

Intense  research  is  also  focused  on  multi-dimensional  finite- volume  upwind  schemes  that  induce 
upwinding  along  a  few  significant  directions.  An  early  effort^  generated  a  grid-independent  upwind 
scheme  based  on  directional  upwinding  along  possible  shock  wave  directions.  This  approach  later 

*  Associate  Professor,  jiannell@utk.edu.  Joe  lannelli,  Dept,  of  Mechanical  Aerospace  Engineering  and  Engi¬ 
neering  Science,  The  University  of  Tennessee,  315  Perkins  Hall,  Knoxville,  TN  37996-2030,  U.S.A. 
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enjoyed  addition  of  local  Riemann  solutions^®  along  several  upwind  directions  including  the  flow 
velocity,  speed  gradient  and  pressure  gradient  directions.  An  alternative  second  order  rotated 
upwind  scheme  used  flux-  diff'erence  splitting  along  two  orthogonal  directions'^  determined  on  the 
basis  of  the  pressure  gradient.  Other  approaches  involved  approximate  multi-dimensional  Riemann 
solvers  and  local  wave  decompositions,  with  wave  modeling.^^"^^  In  these  formulations,  some  wave 
directions  and  strengths  must  still  be  fixed  a-priori  to  generate  a  viable  CFD  algorithm. 

Difficulties  exist  in  these  methods  in  assessing  the  magnitude  of  the  induced  multi-dimensional 
upwind  diffusion.  It  is  also  difficult  to  determine  whether  consistent  upwinding  exists  not  only 
over  the  selected  directions,  but  along  all  flow-field  wave-  propagation  directions.  Additional  data 
filtering  or  upwind-direction  freezing  may  also  be  required  for  convergence  and  essential  mono¬ 
tonicity.  More  fundamentally,  current  multi-dimensional  upwind  schemes  are  recognized^®  to  rest 
upon  much  less  theoretical  support  than  their  one-dimensional  counterparts. 

This  report  expounds  the  multi-dimensional  formulation  of  the  one-dimensional  acoustics  -  con¬ 
vection  upstream  resolution  algorithm.^^  Developed  for  the  2-D  and  3-D  Euler  and  Navier-Stokes 
equations  with  general  equilibrium  equations  of  state,  this  formulation  develops  the  upstream- 
bias  approximation  directly  at  the  differential  equation  level,  before  any  discretization,  and  it 
encompasses,  unifies  and  generalizes  upwind  algorithms,  including  Flux  Vector  Splitting  and  Flux 
Difference  Splitting  developments.  The  formulation  results  in  a  “companion”  characteristics-bias 
system  that  is  associated  with  the  Navier-Stokes  equations  and  rests  on  a  decomposition  of  the 
multidimensional  Euler  Jacobian  into  acoustics  and  convection  components.  Heed  in  particular  that 
no  single  decomposition  of  the  Euler  flux  components  themselves  can  contain  separate  components 
that  respectively  correspond  to  the  physics  of  multi-dimensional  acoustics  and  convection. 

This  formulation  induces  the  upstream  bias  along  all  flow-field  directions  of  wave  propagation 
and  enjoys  a  consistent  theoretical  support  that  rests  upon  the  mathematics  and  physics  of  multi¬ 
dimensional  characteristic  wave  propagation.  The  formulation,  moreover,  reduces  to  a  consistent 
upstream  approximation  of  the  acoustics  equations,  for  vanishing  Mach  number,  which  addresses 
the  challenging  problem  of  calculating  low-Mach-number  flows.  The  formulation,  furthermore, 
supplies  a  stable  and  intrinsically  infinite  directional  upstream-bias  that  induces  minimal  diffusion, 
for  crisp  oblique-  shock  capturing  and  accurate  resolution  of  vortical  structures,  naturally  incorpo¬ 
rates  a  finite  element  discretization  with  implicit  time  integration  because  it  is  straightforward  to 
determine  the  required  Jacobian  matrices. 

A  traditional  centered  discretization  on  arbitrary  grids  of  the  characteristics-bias  system  au¬ 
tomatically  and  directly  generates  a  consistent,  discretely  conservative  and  genuinely  multi  -  di¬ 
mensional  upstream-bias  approximation  of  the  Euler  and  Navier-Stokes  equations.  The  associated 
discrete  multi  -  directional  upstream-bias  remains  independent  of  the  direction  of  the  coordinate 
axes  as  well  as  orientation  of  each  computational-cell  side,  which  obviates  the  need  for  rotated  sten¬ 
cils.  This  approximation  requires  data  only  from  the  computational  cells  shared  by  each  grid  node 
and  also  reduces  to  a  consistent  upstream  approximation  of  the  acoustics  equations,  for  vanishing 
Mach  number,  which  addresses  the  challenging  problem  of  calculating  low-Mach-number  flows. 
Finite  difference,  volume,  or  element  procedures  can  be  used  to  discretize  the  characteristics-bias 
system.  The  algorithm  in  this  paper  has  used  a  finite  element  discretization. 

For  any  Mach  number,  the  characteristics-bias  system  induces  consistent  upwinding  along  all 
wave-propagation  directions  radiating  from  any  flow-field  point.  The  upstream  directions  can 
be  continuously  updated  and  high-rate  convergence  to  machine  zero  achieved  without  additional 
shock-capturing  terms,  data  filtering  or  loss  of  essential  monotonicity. 

The  formulation  induces  an  anisotropic  variable-strength  upstream  bias  that  directly  correlates 
with  the  multi-dimensional  spatial  distribution  of  characteristic  velocities.  The  magnitudes  of  the 
associated  streamwise  and  crossflow  dissipations  remain  different  from  and  independent  of  each 
other,  with  crossflow  dissipation  that  decreases  for  increasing  Mach  number.  In  this  manner  the 
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developed  formulation  will  not  generate  for  increasing  Mach  number  as  much  crosswind  dissipation 
as  induced  by  an  isotropic  or  direction-split  formulation.  The  magnitude  of  the  induced  upstream 
bias,  furthermore,  depends  on  local  solution  smoothness.  Only  at  solution  discontinuities,  e.g. 
a  shock  wave,  is  the  induced  upstreaming  commensurate  with  a  fully  upwind  formulation.  In 
regions  of  solution  continuity,  the  upwind-bias  reduces  to  a  minimum,  which  corresponds  to  minimal 
induced  dissipation. 

The  operation  count  for  this  algorithm  is  comparable  to  that  of  a  simple  flux  vector  splitting 
algorithm.  The  developments  in  this  study  have  employed  basic  four-noded  cells.  To  determine  the 
ultimate  accuracy  of  bi-linear  approximations  of  fluxes  within  four-noded  cells,  for  a  computation¬ 
ally  efficient  implementation,  this  study  employs  no  MUSCL-type  local  extrapolation  of  dependent 
variables. 

This  report  is  organized  in  11  sections  and  7  appendices.  After  the  introductory  remarks  in  Sec¬ 
tion  1,  Section  2  presents  the  governing  equations  and  Sections  3  introduces  the  multi-dimensional 
non-discrete  upstream-bias  formulation.  A  multi-dimensional  characteristics  analysis  for  general 
equations  of  state  is  the  subject  of  Section  4,  followed  in  Section  5  by  the  determination  of  the 
multi-dimensional  acoustics  and  convection  components  within  the  Euler  flux  Jacobian.  Section 
6  presents  the  acoustics-convection  characteristics  Euler  flux  and  Section  7  shows  that  this  for¬ 
mulation  is  genuinely  multi-dimensional.  Section  8  presents  the  finite  element  discretization  of 
the  characteristics-bias  system  and  Section  9  details  the  controller  of  induced  upstream  dissipa¬ 
tion.  Section  10  delineates  a  non-linearly  stable  implicit  Runge-Kutta  algorithm,  with  concluding 
remarks  presented  in  Section  11.  Appendices  A-F  delineate  several  ancillary  mathematical  devel¬ 
opments  and  Appendix  G  contains  a  copy  of  an  article  on  the  1-D  version  of  this  procedure,  as 
published  in  the  International  Journal  for  Numerical  Methods  in  Fluids. 


2  Governing  Equations 


2.1  Navier-Stokes  System 

With  respect  to  an  inertial  Cartesian  reference  frame  and  implied  summation  on  repeated  subscript 
indices,  the  Navier-Stokes  conservation  law  system^^  is 

dt  dxj  dxj 

which  consists  of  the  continuity,  momentum  and  total-energy  equations.  For  3-D  flows,  subscript 
j  varies  in  the  range  1  <  j  <  3;  with  TZ  denoting  the  real-number  field,  the  independent  variable 
{x,t)  =  {xi,X2,xz,t)  in  (1)  varies  in  the  domain  D  =  Q,  x  [to,T],  [to,T]  €  C  1Z^.  The 

dependent- variable  array  q  =  q{x,  t)  and  viscous  and  inviscid  flux  “vector”  components  fj  =  fj  (q) 
and  fj  =  fj{q)  are  then  defined  as 


N 

P 

0 

TUi 

Di 

m2 

>  ,  f]{Q)  =  < 

m3 

^3j 
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fi(g)  = 


-mi  +p 


,  /2(g)  = 


-m2  +p 


,  /3(g)  = 


-m3  -hp 


■(E  +  p) 


'-(E-hp) 


’■(E-hp) 


In  the  array  g,  the  variables  p,  mi,  m2  ,  m3,  and  E,  respectively  denote  static  density  and  volume- 
specific  linear  momentum  components  and  total  energy.  Concerning  the  viscous  and  inviscid  fiuxes, 
the  variables  Tij,  qj  and  p  respectively  indicate  the  deviatoric  stress- tensor  components,  the  Fourier 
heat  conduction  flux  component  and  static  pressure.  The  eulerian  flow  velocity  u,  with  cartesian 
components  Uj,  \  <  j  <  3,  is  then  defined  as  u  =  m/p. 

2.2  Constitutive  Relations 

The  components  Tij  of  the  deviatoric  stress  tensor  are  expressed  as 


Tij  =  p, 


dui  duj 
dxj  dxi 


.  dui  2 

-f  A-5 — Oij  ,  A  —  —-p-{-r)B 


where  p  and  A  respectively  denote  the  first  and  second  coefficient  of  viscosity,  and  tjb  indicates 
bulk  viscosity.  For  monoatomic  gases  Pb  vanishes,  while  for  other  ffuids,  like  air,  Stokes’  hypothesis 
A  =  —  |/i  constitutes  a  reliable  approximation. 

The  components  of  the  Fourier  heat  flux  vector  are  expressed  as 

=  (5) 

where  T  indicates  static  temperature  and  k  denotes  the  coefficient  of  thermal  conductivity. 


2.3  Equations  of  State  and  Speed  of  Sound 

For  any  homogeneous  equilibrium  gas,  pressure  can  be  expressed  as  a  function  of  two  other  ther¬ 
modynamic  variables.^^  They  are  density  p  and  mass-specific  internal  energy  e,  in  this  case,  since 
they  are  readily  available  from  the  Euler  system  (1):  p  directly  from  the  continuity  equation  in  the 
system,  and  e  from  q  as 

+  (6) 

The  pressure  equation  of  state  thus  becomes 


P  =  P{P,  f)  =  P  ^  ("i?  +  m|  -f  mfj 


According  to  this  expression,  the  Jacobian  derivatives  of  p  with  respect  to  q,  for  the  Jacobian  dfj/dq 
of  fj{q),  are  not  all  independent  of  one  another.  The  Jacobian  derivatives  of  (7)  with  respect  to 
mi,  m2,  m3  and  E  in  fact  satisfy  the  constraints 


mi  dp 

p  dE 


m2  dp 
p  dE 


m3  dp 

p  dE 


Acoustics-Convectiop  Upstream  Resolution  Algorithm 


7 


as  obtained  by  expressing  the  derivatives  of  p  with  respect  to  mi,  m2  ,  m3  and  E  in  terms  of  the 
thermodynamic  derivative  ofp  with  respect  to  e,  from  the  first  expression  in  (7).  In  the  following 
sections,  for  simplicity,  the  abridged  notation 


Pp  Qp  >  Pmi  - 


_  ^P 
dmi 


-  ^  -dp  _  dp  ... 

“  dm2  ’  “  dm3  ’  dE 


will  denote  the  Jacobian  derivatives  of  pressure.  In  terms  of  these  derivatives,  the  square  of  the 
speed  of  sound  c  for  general  equations  of  state  can  be  expressed^^  as 


-  4  (™1  +  ™2  +  ">3) 

\  H  r 


This  result,  in  particular,  allows  expressing  the  mass-specific  total  enthalpy  H  as 


_  E  +  p  _  ^  (i+Pe^^)  -Pp 


where  M  =  ||ii||/c  denotes  the  Mach  number. 

The  specific  perfect-gas  expressions  for  (6)  and  (7)  follow  from  the  internal  energy  and  equation 
of  state 

R 

e  =  CvT  = - -T  ,  p  =  pRT  (12) 

7-1 

for  this  type  of  gas,  where  Cy,  T,  R  and  7  respectively  denote  the  constant- volume  specific  heat, 
static  temperature,  gas  constant  and  specific-heat  ratio.  The  elimination  of  T  from  these  two 
expressions,  along  with  (6) ,  leads  to  the  following  familiar  expressions  for  the  equation  of  state  for 
p  in  terms  of  q 

p  =  (7-l)pe  =  (7-1)  -  ^(nil  +  ml+mfj^  (13) 

With  this  equation  of  state,  the  corresponding  Jacobian  derivatives  of  pressure  become 


Pb  =  (7  -  1) 


which  satisfy  (8).  Prom  (10),  the  perfect-gas  speed  of  sound  can  be  expressed  as 


=  (7  -  1) 


E  +  p  +  m2  +  m3 
P  2p2 


directly  in  terms  of  the  dependent  variables. 


3  Non-Discrete  Characteristics-Bias  Approximation 

The  non-discrete,  i.e.  continuum  or  before  discretization,  characteristics-bias  approximation  of  the 
Navier-Stokes  equations  is  a  system  of  equations  that  encompasses  the  Navier-Stokes  equations 
(1)  and  also  induces  an  authentically  multi-dimensional  solution-dependent  upstream  dissipation 
along  not  just  select  directions  but  all  characteristic  directions.  This  upstream  dissipation  originates 
from  a  multidimensional  difierential  parabolic  perturbation  within  the  characteristics-bias  system; 
this  parabolic  perturbation,  in  turn,  emerges  from  a  decomposition  of  the  Jacobian  matrix  of  the 
multi-dimensional  Euler  flux  fj. 
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The  development  of  a  characteristics-bias  approximation  at  the  differential  equation  level, 
prior  to  the  discretization  in  space,  possesses  distinctive  advantages.  For  instance,  the  direc¬ 
tions  of  this  multi-dimensional  upstream  bias  remain  independent  of  any  grid-line  direction,  and 
correlate  with  characteristic  directions;  furthermore  a  straightforward  centered  approximation  of 
the  characteristics-bias  system  automatically  generates  a  genuinely  multi-dimensional  consistent 
upstream-bias  approximation  of  the  original  Navier-Stokes  equations,  without  any  need  for  addi¬ 
tional  numerical  dissipation  terms.  This  formulation,  therefore,  induces  but  a  minimal  amount  of 
upstream  dissipation,  that  amount  that  leads  to  practically  non-oscillatory  accurate  solutions. 

Traditional  discrete  approximations  of  (1)  originate  from  the  integral  statement 


dfjjq) 

dxj 


dil  =  0 


(16) 


which  is  equivalent  to  the  governing  system  (1)  for  arbitrary  subdomains  C  and  arbitrary 
integrable  test  functions  w  with  compact  support  in  For  finite  volume  and  element  formula¬ 
tions,  some  upstream-bias  approximations  can  emerge  from  this  statement  through  numerical  flux 
formulae  for  fj  and  select  choices  for  the  test  functions. 

The  characteristics-bias  formulation  rests  on  the  characteristics-bias  system 


dq  dffjg)  dfVjq)  ^  ^ 

dt  dxj  dxj 

where  /p  corresponds  to  a  characteristics  flux  that  automatically  induces  a  multi-dimensional  and 
infinite  directional  upstream-bias  approximation  for  the  Euler  flux  divergence  The  associated 
characteristic-bias  integral  is  then  defined  as 


dn  =  o 


(18) 


Since  the  characteristics  flux  is  developed  independently  and  before  any  discretization,  a  genuinely 
multi-dimensional  upstream-bias  approximation  for  the  Navier-Stokes  equations  (1)  on  arbitrary 
grids  directly  results  from  a  straightforward  centered  discretization  of  this  integral  statement  on 
any  given  grid.  For  finite  element  formulations,  a  multi-dimensional  upstream-bias  approximation 
of  the  Navier-  Stokes  equations  emerges  from  the  most  basic  Galerkin  approximation  of  (18). 

The  following  sections  show  that  the  divergence  of  the  characteristics  flux  /j^  naturally  derives 
from  an  upstream-bias  integral  average  as 


e=i 


dq 

dxj 


(19) 


In  this  expression,  Aij  corresponds  to  a  matrix  component  of  the  Euler  flux  Jacobian,  such  that  the 
matrix  AijUj  has  uniform-sign  eigenvalues,  where  rij  denotes  the  component  of  a  unit  vector 
along  an  arbitrary  wave-propagation  direction  within  conical  regions  within  the  flow  field.  The 
coefficients  denote  linear-combination  functions,  possibly  depending  upon  q.  The  term  an  in¬ 
dicates  the  direction  cosine  of  a  unit  vector  ai  along  the  principal  wave-propagation  direction 
of  wave  a  convection  or  acoustic  wave  for  instance.  The  positive  V’  in  (22),  0  <  ip  <1,  stands 
for  a  new  “upstream-bias”  controller,  detailed  in  Section  9,  which  controls  the  amount  of  induced 
upstream-bias  dissipation.  This  controller  adjusts  this  dissipation  depending  on  local  solution  non¬ 
smoothness.  In  regions  of  discontinuous  solution  slopes,  this  controller  increases,  to  preserve  numer¬ 
ical  stability  and  capture  shocks  crisply;  in  regions  of  smooth  solutions,  the  controller  approaches 


AcousticS’’Convection  Upstream  Resolution  Algorithm 


9 


its  minimum,  to  reduce  upstream  dissipation.  In  this  manner,  the  characteristics-bias  formulation 
generates  minimal-upstream-dissipation  solutions.  The  positive  e  denotes  a  local  length  scale  to 
ensure  consistency  of  the  characteristics-bias  formulation.  It  is  well  known  that  any  discrete  ap¬ 
proximate  solution  of  the  Navier-Stokes  equations  exactly  solves  not  the  Navier-Stokes  equations, 
but  an  associate  system  of  equations,  also  called  the  modified  equations;  as  a  maximum  grid  spac¬ 
ing  within  the  discrete  approximation  approaches  nought,  the  modified  equations  approach  the 
Navier-Stokes  equations  and  the  discrete  solution  is  expected  to  approach  the  exact  solution  of 
these  equations.  Within  each  cluster  of  nodes  of  a  discretization  grid,  therefore,  e  in  (18)  is  set 
equal  to  a  local  grid  spacing  so  that  the  discrete  solution  of  (18)  is  expected  to  approach  a  solution 
of  the  Navier-Stokes  equations  as  e  approaches  nought. 

The  characteristics  flux  depends  on  physical  characteristic  directions.  Section  4,  therefore, 
details  a  multi-dimensional  characteristics  analysis,  which  clarifies  the  interaction  between  acoustic 
and  convection  waves  for  all  Mach  numbers.  Section  5  determines  the  acoustics  and  convection 
components  within  the  Jacobian  matrices  of  the  Euler  flux.  Section  6  then  presents  the  acoustics- 
convection  form  of  the  characteristics  flux  divergence  (19). 


3.1  Flux  Jacobian  Decomposition  and  Upstream-Bias  Integral  Average 


To  develop  the  flux  /^,  consider  first  the  flux  Jacobian  decomposition  (  FJD  )  into  L  contributions 


dfj  _  .  dq 


(20) 


where  ag  denotes  a  linear-combination  function,  possibly  depending  upon  q,  Agj  corresponds  to  a 
flux-jacobian  matrix  component  such  that  the  matrix  Agjrij  has  uniform-sign  eigenvalues  within  a 
conical  region  spanned  by  Uj,  within  the  flow  domain. 

An  integral  average  of  the  Euler  flux  divergence  ^  as  expressed  through  decomposition  (20) 
becomes 

^  dq 


L  dft=  dO, 

Jn  oxj  Jn  ^  dxj 


(21) 


The  flux  is  therefore  deflned  by  way  of  an  upstream-bias  integral  average  as 


r  df^  r  ^ 

I  dfl=  L^{w  +  -tpSiw)  OLgAg. 

Jn  oxj 


dxj 


(22) 


where  the  rhs  integral  provides  an  upstream  bias  for  each  matrix  component  within  the  FJD  in 

(20). 

The  positive  ijj  in  (22),  0  <  <  1,  stands  for  a  new  “upstream-bias”  controller,  which  auto¬ 

matically  adjusts  the  amount  of  induced  upstream-bias  diffusion,  depending  on  local  solution  non¬ 
smoothness.  The  expression  5gw  denotes  a  directional  variation  of  the  test  function  w  along  the  axis 
of  a  conical  region  within  the  flow  domain.  This  variation  induces  the  appropriate  upstream-bias 
for  the  test  function  w  for  each  component  within  (22).  Depending  on  the  physical  signifi¬ 
cance,  magnitude  and  algebraic  sign  of  the  eigenvalues  of  AgjUj,  the  variation  6gw  can  vanish  or 
become  algebraically  positive  or  negative,  which  corresponds  to  an  upstream  bias  respectively  in 
the  negative  or  positive  sense  of  the  axis  of  each  conical  region. 


3.2  Characteristics-Bias  Flux 


The  directional  variation  6gw  in  (22)  becomes 

-  ^  dw .  dw 
dgw  =  -^dgXi  =  -^aue 


^g^i  —  0/ig£ 


(23) 
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where  e  denotes  a  local  positive  length  scale  and  an  indicates  the  z**'  direction  cosine  of  a  unit 
vector  a£  along  the  principal  wave-propagation  direction  of  wave  ‘T’ . 

With  these  specifications,  the  upstream-bias  integral  average  (22)  becomes 


n  ox 


^dfi 


(24) 


Considering  that  w  has  compact  support  in  fl,  it  vanishes  on  the  boundary  dQ  of  fl.  As  a  result, 
integrating  (24)  by  parts  generates 


e^Y^anaeAtj-^ 

.  e=i  ^^3 ) 


dQ,  =  0 


(25) 


which  contains  no  boundary  integrals.  Since  this  integral  must  vanish  for  arbitrary  test  functions  w 
and  domains  !!1,  its  integrand  must  identically  equal  zero,  which  generates  the  following  expression 
for  the  divergence  of  the  characteristics  flux 


dxi 


dx 


*  \  ^=1  ^^3 


(26) 


This  expression  exhibits  an  upstream-bias  artificial  diffusion,  in  the  form  of  a  second-order  differ¬ 
ential  expression  with  associated  upstream-bias  matrix 


A  = 


(27) 


where  rii  indicates  the  direction  cosine  of  a  unit  vector  n  along  an  arbitrary  wave-propagation 
direction.  Appendix  D  shows  that  for  physical  consistency  of  the  upstream  bias  in  (22),  (19)  and 
associated  mathematical  stability  of  the  corresponding  second-order  differential  expression,  all  the 
eigenvalues  of  this  upstream-bias  matrix  must  be  positive  at  every  flow-field  point  and  for  any  wave- 
propagation  direction  n.  This  requirement  implies  a  correct  upstream  along  all  directions  radiating 
from  any  flow  field  point  and  thus  becomes  a  fundamental  upstream-bias  stability  condition. 

Prom  (19),  an  expression  for  the  component  of  the  characteristics  flux  is 


f?  =  fi-ei£j^aaaiAii^  (28) 

As  a  multi-dimensional  expression,  each  cartesian  component  //"  also  depends  on  the  derivatives 
of  the  solution  q  along  all  cartesian  coordinate  directions.  The  continuum  expression  (19),  or  (28), 
thus  constitutes  a  non-discrete  generalization  of  the  various  numerical-flux  formulae  employed 
in  several  CFD  upwind  schemes.  It  encompasses,  generalizes,  and  unifies  flux-vector  and  flux- 
difference  schemes  as  shown  by  the  following  representative  1-D  examples. 


3.2.1  van  Leer’s  Formulation  and  Flux  Vector  Splitting 

Consider  the  van  Leer’s  formulation  as  a  representative  Flux  Vector  Splitting  (  FVS  ).  In  this 
formulation,  the  inviscid  flux  /  is  “split”  as 


/  =  + /«- 


(29) 
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where  the  Jacobian  matrices  of  and  respectively  possess  non-negative  and  non-positive 
eigenvalues. 

The  FJD  expression  (20)  encompasses  (29)  with  L  =  2  as 

L  QfVL+  QfVL- 

+  .  “1  =  1  .  “2  =  1  (30) 

The  corresponding  characteristics-bias  flux  divergence  for  van  Leer’s  FVS  accrues  from  (19) 
with  V’  =  1,  oi  =  1,  02  =  -1  as 

df^  _df  d  (  ( \dq\  df  d  (  (dl^_d^\\ 
dx  dx  5a;  y  y  dq  dq  J  dx  J  dx  dx 

which  generalizes  in  the  continuum  the  traditional  numerical  flux  formulae  for  FVS  constructions. 
The  associated  upstream  matrix  A  is 


QfVL+  qjvl- 

dq  dq 


(32) 


The  upstream-bias  stability  condition,  however,  is  not  automatically  satisfied,  even  though  each  of 
the  two  matrices  ^  and  ^  has  positive  eigenvalues.  This  stability  condition  is  not 

unconditionally  satisfied  because  the  sum  of  two  positive-eigenvalue  matrices  does  not  necessarily 
yield  a  matrix  with  positive  eigenvalues.  As  an  example  consider  the  following  matrix  sum  of  two 
positive-eigenvalue  matrices 


where  cr  is  a  real  number.  One  of  the  eigenvalues  of  this  matrix  sum  is  negative  for  a  >  4.  For 
instance,  for  cr  =  7  the  eigenvalues  are  -1-9  and  —1. 

Most  likely,  however,  (32)  satisfies  the  upstream-bias  stability  condition  for  most  of  the  flow 
conditions  considered  in  the  technical  literature,  in  view  of  the  stable  results  reported.  For  sub¬ 
sonic  flows,  each  of  the  two  flux  vector  components  in  (29)  remains  unrelated  to  the  physics  of 
acoustics  or  convection.  On  the  other  hand,  (29)  is  computationally  advantageous,  for  it  calls  for 
the  discretization  of  simple  flux-vector  components. 


3.2.2  Roe’s  Formulation  and  Flux  Difference  Splitting 

Consider  next  Roe’s  formulation  as  a  representative  Flux  Difference  Splitting  (  FDS  )  development. 
In  this  formulation,  the  inviscid  flux  Jacobian  of  /  is  “split”  as 

^  -h  XA- (34) 

where  X  and  A  =  A"'‘4-A“  denote  the  right  eigenvector  matrix  and  eigenvalue  diagonal  matrix  of  the 
Jacobian,  all  evaluated  at  special  average  values  of  q,  with  A"*"  and  A~  respectively  containing  non¬ 
negative  and  non-positive  eigenvalues.  The  matrices  at  the  rhs  of  (34),  therefore,  will  respectively 
possess  non-negative  and  non-positive  eigenvalues. 

The  FJD  expression  (20)  encompasses  (34)  with  L  =  2  as 

L 

Y,  =  XA+X-^  +  XA-X-1 


,  Q-l  =  1  ,  Q!2  =  1 


(35) 
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The  corresponding  characteristics-bias  divergence  for  Roe’s  formulation  accrues  from  (19)  with 
=  1,  oi  =  1,  02  =  — 1  as 

d 

dx  dx 


I  (.  (XA-Jf-  -  XA-X-‘)  I)  =  g  -  «  (.X  (A-  -  A- )  X->|)  (36) 


which  generalizes  in  the  continuum  the  traditional  numerical  flux  formulae  for  FDS  constructions. 
The  associated  uspstream  matrix  A  is 


=  X  ( A+  -  A- )  X 


-1 


(37) 


which  has  non-negative  eigenvalues  and  therefore  automatically  satisfles  the  upstream-bias  stability 
condition  for  any  flow  state  for  which  no  eigenvalue  vanishes.  The  discretization  of  (36)  calls  for 
more  computational  operations  than  (31),  while  each  of  the  two  rhs  components  in  (34)  lumps  into 
one  matrix  the  matrices  representative  of  the  distinct  acoustics  and  convection  wave  propagation 
mechanisms.  On  the  other  hand,  numerous  numerical  results  bear  out  the  accuracy  of  an  FDS 
formulation. 


4  Characteristics  Analysis 

Within  a  flow  field,  acoustic  and  convection  waves  propagate  in  infinitely  many  directions,  and 
along  each  direction  the  associated  propagation  velocity  also  depends  on  the  Mach  number.  This 
section  introduces  a  new  intrinsically  multi-dimensional  characteristics  analysis  based  on  a  non¬ 
linear  wave-like  form  of  the  solutions  of  the  Euler  equations.  This  analysis  leads  to  the  spatial 
distribution  of  multi-dimensional  propagation  velocities  and  shows  that  among  all  propagation 
directions  the  streamline  and  crossflow  directions  are  principal  propagation  directions.  This  line  of 
inquiry  yields  specific  conditions  for  a  physically  coherent  upstream  bias  formulation  that  remains 
consistent  with  multi-dimensional  acoustic  and  convection  wave  propagation.  Several  details  are 
presented  in  Appendix  B,  which  also  delineates  a  characteristic-surface  analysis  and  indicates  that 
the  streamline  and  crossflow  directions  are  principal  directions  of  wave  propagation  in  the  time- 
space  continuum. 

4.1  Non-linear  Wave-Like  Solutions 

With  implied  summation  on  repeated  indices,  the  non-linear  wave-like  form  of  solutions  q  of  the 
Euler  equations  is  expressed  as 

q  =  q{r}i)  ,  t]\  =  x  n  —  X{q)t  =  XjUj  —  X{q)t  (38) 

where  n  denotes  a  space-domain  propagation-direction  unit  vector,  independent  of  and 

A  =  X{q)  indicates  a  wave-propagation  velocity  component  along  the  n  direction.  Appendix  A 
delineates  many  important  properties  of  non-linear  wave-like  solutions. 

4.2  Characteristic  Velocity  Components 

The  solution-dependent  velocity  component  A  =  A(g)  is  determined  by  enforcing  the  condition  that 
q  as  expressed  in  (38)  satisfies  the  Euler  equations.  Appendix  A  shows  that  this  condition  yields 
the  eigenvalue  problem 


(39) 
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For  non-trivial  solutions  hence  q  =  q{T]i),  the  characteristic  velocity  components  A  are  thus 
the  eigenvalues  of  the  linear  combination  of  flux  vector  jacobians 


dfjjq) 

dq 


rij  = 


(40) 


Til 

1  rL2 

1  n3 

,  0  \ 

-ujuj-nj  -f  , 

uini  +  ujuj  ni 

.  wi7i2+p^2^1 

,  Wins  +  p^g  ni 

,  P^jTll  \ 

-\-ppn2  , 

U2ni  +  712 

,  U2T12  +p^2  ”2 

.  U2n3  -1-  p^g  712 

-usUjTi j  -\-p^ns 

usni  +p^^  713 

,  +pTn3»^3 

>  Pb^s  I 

\  {Pp  -  Ji) 

Htii  +ujnjp^^ 

,  Hn2  + 

,  Hns^Ujnjp^^ 

•  “i"j  (l+PE)  / 

For  general  equations  of  state,  these  eigenvalues  have  been  exactly  determined  in  closed  form  as 


UjUj 


,  Xi%  =  ujnj  ujuj) ) 


(41) 


where  superscript  ds  signifies  dimensional  Euler  eigenvalues.  Of  interest,  eigenvalues  directly 
incorporate  a  sound  speed  expression  that  coincides  with  the  isentropic  partial  derivative  of  pressure 
(10).  Through  (10),  therefore,  these  equilibrium-gas  eigenvalues  become 

Xt%  =  Ujnj  ,  A3^  =  Uj-njic  (42) 

which  have  the  same  familiar  form  as  the  perfect-gas  eigenvalues.  The  non-dimensional  form  of 
(42)  follows  firom  division  by  c,  which  supplies  the  Mach-number  dependent  expressions 

Xq  2  =  VjUjM  ,  A3^4  =  vjrijM  ±  1  (43) 

where  vi,  V2  and  vz  denote  the  components  of  a  unit  vector  v  in  the  velocity  u  direction. 

As  the  inner  product  of  the  two  unit  vectors  n  and  u,  the  contraction  vjrij  supplies  the  cosine 
of  the  angle  {9  —  0^,)  between  n  and  u,  where  6  and  respectively  denotes  the  angle  between  n 

and  the  xi  axis  and  the  angle  between  v  and  the  xi  axis  on  any  plane  11  that  contains  the  velocity 

vector.  The  eigenvalues  (43)  thus  become 

Ao_2  =  cos(0  — 0t,)Af  5  A3  4  =  cos  (0  —  ^v)  Af  ±  1  (44) 

These  expressions,  in  particular,  imply  that  the  Euler  eigenvalues  achieve  their  extrema  for  0 
hence  when  n  points  in  the  streamline  direction,  whereas  for  n  pointing  in  the  crossflow  direction, 
hence  6  =  90°  +  6v,  these  eigenvalues  no  longer  depend  upon  M. 

The  convection  eigenvalues  Ag  2  vanish  when  cos  {6  —  0^,)  =  0,  hence  for  n  perpendicular  to  the 
streamline  direction,  or,  equivalently,  pointing  in  the  crossflow  direction.  Since  ||  cos  {9  —  9^,)  ||  <  1, 
the  acoustic-convection  eigenvalues  A3  4  can  only  vanish  for  Af  >  1,  hence  for  supersonic  flows.  For 
these  flows,  A3  4  =  0  for 


T  cos  {9  -9^)  =  ±  sin  ((0  -  90°)  -9^)  =  ^  (45) 

hence  for  n  pependicular  to  the  Mach  lines,  because  ±((0  —  90°)  —  9-o)  corresponds  to  the  angle 
between  a  Mach  line  and  the  streamline,  from  the  well  known^®  second  expression  in  (45).  The 
lines  that  are  perpendicular  to  the  Mach  lines  will  be  called  “conjugate”  lines. 

The  lines  that  are  perpendicular  to  n  for  vanishing  eigenvalues  Ag  2  and  A3  4  thus  respectively 
become  the  streamline  and  Mach  lines.  This  result  is  not  coincidental,  for  vanishing  eigenvalues 
Ai  4  correspond  to  wave-like  solutions  of  the  steady  Euler  equations,  for  which  the  streamline  and 
Mach  lines  are  characteristic-wave  propagation  lines. 
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4.3  Polar  Variation  of  Characteristic  Speeds 

Figure  1  presents  the  polar  variation  of  the  absolute  values  of  eigenvalues  (43)  for  subsonic,  sonic 
and  supersonic  Mach  numbers,  in  a  neighborhood  of  a  flow  field  point  P  on  any  plane  11  that 
contains  the  velocity  vector  at  P. 


M  =  0.05 

r.  100°  PO®  80° 

^  110°  70° 


M  =  0.5 

^  100^  90^  80^ 


70° 
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Figure  1:  Polar  Variation  of  Subsonic  Wave  Speeds 


M  =  1.5 


M  =  2.0 


Figure  2:  Polar  Variation  of  Supersonic  Wave  Speeds 
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These  variations  are  obtained  for  a  variable  unit  vector  n  =  (cos  d,  sin  6)  and  fixed  unit  vector 
V,  in  this  representative  case  inclined  by  +30°  with  respect  to  a  reference  Xi  axis  on  11.  A  collective 
inspection  of  all  these  diagrams  reveals  three  shared  features  for  all  Mach  numbers.  The  maxi¬ 
mum  characteristic  speeds  occur  in  the  velocity  direction,  i.e.  along  a  streamline,  as  noted  before. 
Secondly,  all  the  characteristic  speeds  are  symmetrically  distributed  about  the  streamline  direc¬ 
tion.  Thirdly,  the  eigenvalue  pairs  (||Af  ||,  HAf  ||)  and  (||Af  ||,  IIA4  |1)  remain  mirror  skew-symmetric 
with  respect  to  the  crossflow  direction,  in  the  sense  that  the  curves  representative  of  HA^  ||  and 
II A4  II  become  the  respective  mirror  images  of  the  variations  of  ||Ai  ||  and  IIA3  ||  with  reference  to 
this  direction.  The  streamline  and  crossflow  directions,  therefore,  become  two  fundamental  wave- 
propagation  axes. 

For  vanishing  Mach  numbers,  the  acoustic-convection  propagation  curves  in  the  figure  approach 
two  circumferences.  The  distribution  of  propagation  speeds  in  this  case  is  therefore  isotropic,  which 
corresponds  to  the  direction-invariant  propagation  of  acoustic  waves.  As  the  Mach  number  increases 
from  zero,  the  curves  in  the  figure  progressively  become  circular  asymmetric,  which  corresponds  to 
anisotropic  wave  propagation.  For  M  =  0.5  this  anisotropy  is  already  evident  and  becomes  more 
pronounced  for  higher  Mach  numbers.  The  non-dimensional  characteristic  speeds  then  approach  1 
in  the  region  about  the  crossflow  direction,  which  corresponds  to  essentially  acoustic  propagation. 

For  all  Mach  numbers,  the  convection  eigenvalues  Ao,2  change  sign  when  the  n-direction  shifts 
from  an  upstream  to  a  downstream  axis  with  respect  to  u.  For  this  reason,  the  associated  curves 
cross  the  polar  origin.  Pure  convective  propagation,  therefore,  remains  mono-axial,  firom  upstream 
to  downstream  of  P,  and  the  axis  of  this  type  of  wave  propagation  is  the  streamline. 

For  subsonic  Mach  numbers  the  eicoustic-convection  eigenvalues  Af  and  A4  respectively  remain 
positive  and  negative  for  all  directions.  For  this  reason  the  associated  curves  contain  the  polar 
origin.  For  subsonic  flows,  therefore,  acoustic-convection  waves  propagate  bi-modally,  from  both 
upstream  and  downstream  toward  and  away  firom  point  P,  along  all  directions  radiating  fi:om  P. 

Beginning  at  the  sonic  state,  this  pattern  drastically  changes  for  supersonic  Mach  numbers.  In 
this  case  both  A3  and  A4  change  algebraic  sign  when  the  n  shifts  from  upstream  to  downstream  of 
P  along  a  streamline.  For  this  reason,  the  associated  curves  cross  the  polar  origin.  For  supersonic 
flows,  therefore,  acoustic-convection  wave  propagation  becomes  mono-axial  along  a  streamline, 
firom  upstream  to  downstream  of  P. 

4.4  Regions  of  Supersonic  and  Subsonic  Wave  Propagation 

The  streamline  and  crossflow  directions  are  principal  directions  of  gas  dynamic  acoustic  and  con¬ 
vection  wave  propagation.  These  directions  become  the  axes  of  distinct  wave-propagation  regions, 
which  in  this  report  are  named  the  streamline  and  crossflow  wave  propagation  regions. 

The  supersonic  mono-axial  wave  propagation  pattern  does  not  extend  to  the  entire  plane  H, 
but  remains  confined  within  two  disconnected  wedge  regions  about  the  streamline.  About  the 
crossflow  direction,  in  fact,  there  exist  two  other  wedge  regions  within  which  acoustic-convection 
wave  propagation  remains  bi-modal,  from  both  upstream  and  downstream  toward  and  away  from 
point  P,  even  for  supersonic  flows.  These  regions  are  determined  by  finding  the  lines  where  the 
Euler  eigenvalues  vanish.  These  lines  are  the  crossflow  and  conjugate  lines,  as  shown  for  subsonic 
and  supersonic  flows  in  Figure  3,  along  with  the  bi-modal  and  mono-axial  propagation  regions,  on 
plane  11,  and  the  domain  of  dependence  and  region  of  influence  of  point  P. 


16 


Jop.  lannalli 


a)  M  =  0.5  b)  M  =  2.0 


Figure  3:  Wave  Velocity  Distribution:  a)  Subsonic  Flows,  b)  Supersonic  Flows 

As  Figure  3-a  shows,  subsonic  acoustic-convection  waves  propagate  bi-  modally  over  the  entire 
plane,  hence  the  Euler  eigenvalues  do  not  all  have  the  same  algebraic  sign.  For  supersonic  flows. 
Figure  3-b  shows  the  streamline  regions  of  mono-axial  wave  propagation,  where  the  Euler  eigen¬ 
values  all  have  the  same  algebraic  sign.  The  conjugate  lines  mark  the  boundary  with  the  crossflow 
regions  of  bi-modal  wave  propagation,  where  the  Euler  eigenvalues  do  not  all  have  the  same  al¬ 
gebraic  sign,  like  the  subsonic-flow  situation.  In  particular,  Ao,2  <  0  and  Ao,2  >  0  respectively 
upstream  and  downstream  of  the  crossflow  direction. 

Since  the  conjugate  lines  remain  perpendicular  to  the  Mach  lines,  it  becomes  easy  to  describe  the 
growth  of  the  mono-axial  propagation  regions  as  the  Mach  number  increases.  For  M  <  1,  no  mono- 
axial  propagation  region  exists.  For  M  =  1,  the  Mach  lines  coincide  with  the  crossflow  direction,  the 
conjugate  lines  coincide  with  the  streamline  direction,  hence  the  mono-axial  propagation  region 
consists  of  the  streamline  only.  As  the  Mach  number  increases,  the  Mach  lines  approach  the 
streamline  and  the  conjugate  lines  approach  the  crossflow  direction.  The  mono-axial  propagation 
regions  thus  grow  as  circular  sectors  of  increasing  angle,  and  the  bi-modal  propagation  regions 
simultaneously  shrink.  Only  for  a  theoretically  infinite  M  will  the  mono-axial  propagation  region 
spread  throughout  the  entire  plane  H.  The  presence,  perhaps  unexpected,  of  these  bi-modal- 
propagation  regions  even  for  supersonic  flows  is  readily  reconciled  with  the  existence  of  the  domain 
of  dependence  and  region  of  influence  for  P  by  the  observation  that  the  bi-modal-propagation 
regions  simply  correspond  to  spaces  within  the  domains  of  dependence  and  ranges  of  influence  of 
other  flow  field  points  besides  P. 

4.5  Physical  Multi-dimensional  Upstream  Bias 

Since  for  all  Mach  numbers  there  simultaneously  exist  regions  of  mono-axial  and  bi-modal  prop¬ 
agation,  a  physically  consistent,  intrinsically  multi-dimensional,  and  infinite  directional  upstream 
formulation  for  the  Euler  and  Navier-Stokes  equations  has  to  provide  an  upstream  approxima¬ 
tion  suitable  for  supersonic  flows  within  the  mono-axial  region,  but  consistent  with  subsonic  wave 
propagation  within  the  bi-modal  region.  This  formulation  has  to  involve  a  streamline  upstream  ap- 
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proximation  that  remains  bi-modal  for  subsonic  flows  and  then  becomes  mono-axial  for  supersonic 
flows.  For  all  Mach  numbers,  a  physical  upstream  approximation  must  then  induce  a  bi-modal 
upstream  bias  along  the  crossflow  direction. 

This  upstream  approximation  has  to  introduce  an  upstream  bias  along  all  directions  radiating 
from  each  flow  fleld  point.  The  bias  in  this  approximation  must  change  with  varying  direction  and 
correlate  with  the  directional  distribution  of  the  characteristic  propagation  speeds.  Furthermore, 
the  directional  distribution  of  the  upstream  bias  must  remain  symmetrical  liVe  the  propagation 
speeds  with  respect  to  both  the  crossflow  and  streamline  directions  on  any  plane  11  that  contains 
the  velocity  vector. 

The  formulation  presented  in  this  report  first  identifies  the  genuine  streamline  and  crossflow 
convection  and  acoustics  components  within  the  flux  Jacobian  matrices.  The  formulation  then 
establishes  a  physically  consistent  upstream  approximation  for  each  of  these  components,  along  all 
wave-propagation  directions,  with  an  upstream-bias  magnitude  that  correlates  with  the  directional 
distribution  of  the  characteristic  propagation  speeds.  For  increasing  Mach  numbers,  from  zero  to 
arbitrarily  high,  the  approximation  gradually  changes  from  bi-modal  in  subsonic  flows  to  mono- 
axial  in  supersonic  flows,  within  the  domain  of  dependence  and  range  of  influence  of  each  flow- 
field  point,  while  remaining  bi-modal  for  all  Mach  numbers  within  an  appropriate  region  about  the 
crossflow  direction. 


5  Acoustics- Convection  Flux  Jacobian  Components 

This  section  determines  the  flux  Jacobian  components  that  correspond  to  multidimensional  con¬ 
vection  and  acoustics.  The  characteristic  analysis  of  the  acoustics  components  shows  that  the 
three-dimensional  isotropic  acoustic  propagation  naturally  decomposes  into  three  directional  prop¬ 
agations,  the  first  along  the  streamline  direction  and  the  remaining  two  along  any  two  mutually 
perpendicular  crossflow  directions;  these  two  cross  flow  propagation  ,  in  particular  ,  combine  into 
a  two-dimensional  isotropic  acoustic  propagation  on  any  plane  perpendicular  to  the  velocity  direc¬ 
tion. 


5.1  Convection  and  Pressure-Gradient  Components 

The  flux  divergence  can  be  decomposed  into  convection  and  pressure-gradient  components  as 


dxj 


dxj  dxj 


(46) 


where  fj  and  /J  respectively  denote  the  convection  and  pressure  flux  components,  defined  as 
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For  supersonic  flows,  the  Euler  eigenvalues  (42),  associated  with  all  have  the  same  algebraic 
sign  within  the  streamline  region  and  the  entire  flux  divergence  can  be  upstream  approximated 
along  the  streamline  principal  direction,  within  this  region.  For  subsonic  flows  these  eigenvalues 
have  mixed  algebraic  sign  and  an  upstream  approximation  for  the  flux  divergence  along  one  sin¬ 
gle  direction  remains  inconsistent  with  the  two-way  propagation  of  acoustic  waves.  Without  the 
pressure  gradient  in  the  momentum  equation,  however,  the  corresponding  flux-jacobian  eigenvalues 
all  have  the  same  algebraic  sign®  within  the  streamline  region  and  the  resulting  convection  flux 
divergence  can  then  be  upstream  approximated  along  one  single  direction.  The  flux  divergence  can 
thus  be  decomposed  as  the  linear  combination 


djl.djl 

dxj  ^  dxj 
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dfl 


0<(3<1 


(48) 


where  the  positive  pressure-gradient  partition  function  /3  can  be  chosen  in  such  a  way  that  all  the 
eigenvalues  associated  with  each  of  the  two  components  between  brackets  in  (48)  keep  the  same 
algebraic  sign  within  the  streamline  region,  for  all  Mach  numbers.  In  this  manner,  these  entire 
components  can  be  upstream  approximated  along  the  streamline.  This  choice  for  /3  is  possible 
because  the  eigenvalues  of  a  matrix  are  continuous  functions  of  the  matrix  entries^®  and  hence 
all  the  eigenvalues  for  the  components  in  (48)  will  continuously  depend  upon  (3.  The  function  P 
will  gradually  increase  toward  1  for  increasing  Mach  number,  so  that  an  upstream  approximation 
for  the  components  in  (48)  smoothly  approaches  and  then  becomes  an  upstream  approximation 
for  the  entire  along  the  streamline,  which  also  ensures  consistency  with  1-D  formulations. 
Decomposition  (48)  is  thus  used  for  an  upstream  approximation  of  the  flux  divergence  within  the 
streamline  region,  for  subsonic  and  supersonic  flows. 

Decomposition  (48),  however,  is  insuflicient  for  an  accurate  multi-dimensional  upstream  mod¬ 
eling  of  acoustic  waves,  for  a  large  Mach  number  range  within  the  crossflow  region  and  for  low 
and  vanishing  Mach  numbers  within  the  streamline  region.  Concerning  the  cross-flow  region,  the 
eigenvalues  associated  with  (48)  become  the  Euler  eigenvalues  (42)  for  /?  =  1,  hence  for  supersonic 
flows.  As  indicated  in  Section  4,  these  eigenvalues  do  not  all  have  the  same  algebraic  sign  within 
the  cross-flow  wave  propagation  region.  A  mono-axial  upstream-bias  approximation  of  (48)  within 
the  cross-flow  region,  therefore,  remains  inconsistent  with  the  existing  two-way  wave  propagation 
in  this  region.  Concerning  the  streamline  wave  propagation  region,  for  a  Mach  number  that  ap¬ 
proaches  zero,  the  Euler  eigenvalues  (42)  can  all  keep  the  same  algebraic  sign  only  if  the  sound  speed 
contribution  vanishes,  which  corresponds  to  a  vanishing  pressure  gradient  contribution  and  hence 
P  approaching  zero.®  But  for  P  approaching  zero,  the  eigenvalues  associated  with  the  components 
in  (48)  approach  the  eigenvalues  of  the  jacobians 
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Using  the  pressure-derivative  identity  (8)  the  eigenvalues  of  these  jacobians  respectively 


are 


•^0,3  —  >  -^4  —  +Pb) 


(51) 


and 

^0,3  =  0  .  K  =  -'^jnjPE  (52) 

which  certainly  all  keep  the  same  algebraic  sign  for  any  wave  propagation  direction,  but  for  van¬ 
ishing  Mach  number  remain  far  less  than  the  dominant  speed  of  sound  c.  For  low  Mach  numbers, 
therefore,  an  upstream  approximation  for  the  components  in  (48)  would  inaccurately  model  the 
physics  of  acoustics.  This  difficulty  is  resolved  by  further  decomposing  the  pressure  gradient  in  (48) 
in  terms  of  genuine  streamline  and  crossflow  acoustic  components,  for  accurate  upstream  modeling 
of  acoustic  waves  throughout  the  flow  field. 


5.2  Acoustic  Equations 

Let  {vi,V2,V3)  denote  the  components  of  a  unit  vector  v  locally  parallel  to  the  velocity  u.  Upon 
expressing  the  velocity  components  {ui, 112,113)  in  terms  of  the  Mach  number  M  as 


(«1,U2,U3)  =  cM{vi,V2,V3)  (53) 

and  using  the  pressure  derivative  identities  (8),  the  Euler  system  becomes 
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where  Sj  and  Cj  respectively  denote  Kronecher’s  delta  and  a  completion  matrix.  The  form  of  this 
matrix  is  Cj  = 
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for  which  Aj°  will  denote  the  acoustics  matrix  multiplying  the  gradient  of  q  in  (56).  The  dependent 
variables  in  these  equations  obviously  correspond  to  those  in  a  flow  fleld  that  originates  from  slight 
perturbations  to  an  otherwise  quiescent  field. 

Heed,  in  particular,  that  in  this  case  the  energy  equation  toward  steady  state  is  no  longer 
linearly  independent  from  the  continuity  equation.  Moreover,  the  momentum  equations  become 
linearly  dependent  upon  each  other  for  steady  wave-like  solutions  q  =  q{fji)  =  q{xjnj).  For  these 
solutions,  in  fact,  the  acoustics  momentum  equations  become 


dp  BE 


l<i<3 


(57) 


This  linear-dependence  phenomenon  directly  explains  the  widely  reported  convergence  difficulties 
towards  steady  state  experienced  in  the  CFD  simulation  of  incompressible,  i.e,  low-Mach-number, 
flows  with  a  compressible  flow  formulation.  The  infinite-directional  algorithm  detailed  in  this  report 
will  incorporate  a  convenient  upstream-bias  expression  to  eliminate  within  the  discrete  equations 
this  steady-state  linear-dependence  phenomenon. 

As  Appendix  A  shows,  the  propagation  velocity  components  A  for  a  wave-like  solution  (38)  of 
this  system  correspond  to  the  eigenvalues  of  the  matrix  which  are 
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and  remain  independent  from  the  wave-propagation  direction  unit  vector  n.  Both  c  and  the 
eigenvalues  A3‘4  in  (58)  correspond  to  the  zero- Mach- number  isentropic  speed  of  sound  (10).  With 
Ag®2  =  0,  the  propagation  of  acoustic  waves  governed  by  this  system,  therefore,  corresponds  to  an 
isentropic  process  with  negligible  flow  kinetic  energy. 

Equations  (56)  contain  the  important  partial-derivative  relations 
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which  result  from  the  momentum  equations  by  expressing  p^  therein  using  (11)  as 
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For  an  isentropic  flow  p  =  p(p),  hence  the  first  rhs  term  in  (60)  equals  The  second  rhs  term 
must  consequently  vanish,  which  returns  results  (59). 

On  the  basis  of  (59)  the  acoustics  equations  become 


dp  /dE  E  Ap  dp  \ 

dxi  \5xj  p  dxi) 


(60) 


'  dp  _  druj 

^  dt  dxj 

I  drrii  _  2 

,  dt  ^  dxi 


(61) 


For  wave-like  fields,  these  equations  yield  the  fundamental  result 
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as  obtained  from  expressions  (137)  in  Appendix  A.  Specifically,  the  first  of  (137)  for  mj  and  the 
second  of  (137)  for  p  are  inserted  into  the  momentum  equations  in  (61)  to  substitute  for 
and  The  partial  derivative  ^  in  the  resulting  momentum  equations  is  then  replaced  by  the 
continuity  equation  in  (61),  which  leads  to  (62).  These  expressions  then  simplify  to 

dmi  dmi  dm2  dmz 

-5 — ”2  +  -X — ”3  =  -X — ni  +  — ni 

0x2  oxz  0x2  oxz 


dm2 

dxi 


ni  + 


dmi 

dxi 


n2  + 


dxz 


dmz  dmz 

■X — ni  +  - — n2  = 
dxi  dx2 


dmi  dm2 

-X — nz  +  -X — nz 
dxi  dx2 


(63) 


Results  (59),  (62),  and  (63)  are  used  in  Section  5.3.2  to  simplify  an  upstream-bias  formulation  for 
the  acoustic  components  within  the  Euler  flux  Jacobian. 


5.3  Acoustic  Components 

For  arbitrary  Mach  numbers  and  corresponding  dependent  variables  p,  mi,  m2,  m3,  and  E,  the 
Euler  flux  jacobians  in  (40)  can  be  decomposed  as 


^_?f[  . 

dq 


dq  ^  dq 


-^+A^+Ar 


(64) 


where,  with  reference  to  the  acoustics  equations  (56),  the  matrices  A^  and  are  defined  as 
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Pe 


0  \ 
Pe4 
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0 


/ 

0\ 
0 
0 
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^4 


0 


(65) 


(66) 


Heed,  in  particular  that  no  flux  component  of  fj{q)  exists,  of  which  the  Jacobian  equals  A“.  The 
eigenvalues  of  the  matrix  Aj^Uj  have  been  determined  in  closed  form  as 


\nc  _  fi 
•^0,3  —  w  » 


A”'  =  —cMp^VjUj 


(67) 


which  become  infinitesimal  for  vanishing  M.  The  matrix  can  be  termed  a  “non-linear  coupling” 
matrix,  for  it  completes  the  non-linear  coupling  between  convection  and  acoustics  within  (64)  so 
that  the  two  Euler  eigenvalues  in  (42)  do  correspond  to  the  sum  of  convection  and  acoustic 
speeds.  Since  the  matrix  A®  will  be  used  in  the  upstream-bias  formulation  for  small  Mach  numbers 
only  and  considering  that  the  eigenvalues  in  (67)  vanish  for  these  Mach  numbers  and  also  for 
n  pointing  in  the  crossflow  direction,  for  which  VjUj  =  0,  no  need  exists  to  involve  A"*^  in  the 
upstream-bias  approximation  of  the  flux  Jacobian  (40). 
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The  eigenvalues  of  AjUj  have  been  exactly  determined  in  closed  form  as 

AS.2  =  0  ,  A3“4  =  ±c  (68) 

and  remain  independent  of  the  propagation  vector  n,  which  signifies  isotropic  propagation.  The 
matrix  Aj,  therefore,  can  be  termed  the  “acoustics”  matrix,  for  its  eigenvalues,  unlike  (51)-(52) 
approach  the  speed  of  sound  c  for  decreasing  Mach  number.  This  matrix,  therefore,  can  be  used 
for  an  upstream-bias  approximation  of  the  Euler  equations  in  the  low  Mach-number  regime,  within 
the  streamline  region,  and  for  any  Mach  number,  within  the  crossflow  region. 


5.3.1  Streamline  and  Crossflow  Components 

For  any  three  mutually  orthogonal  unit  vectors  a  =  {ai,  02,03),  =  {o^^ ,02^ ,0^^^),  and  — 

(oi^ ,  op ,  0^^)  within  a  3-D  flow,  along  with  implied  summation  on  repeated  indices,  the  acoustics 
component  within  the  Euler  flux  divergence  can  be  expressed  as 


(69) 


For  a  parallel  to  w,  this  expression  corresponds  to  a  decomposition  of  the  Euler  acoustics  component 
into  one  streamline  and  two  crossflow  acoustics  components.  For  wave-like  solutions  (38),  three 
eigenvalues  of  each  component  vanish;  the  remaining  eigenvalues  of  these  three  separate  components 
have  been  determined  as 


"3,4 


db  cajUj 


AS  =  ^  cap  rij 


AS  =  ^copuj 


(70) 


The  two  non-vanishing  eigenvalues  associated  with  the  entire  acoustics  component  at  the  Ihs  of 
(69),  but  as  expressed  as  the  rhs  combination  of  streamline  and  crossflow  components  have  then 
been  determined  as 


^3,4 


=  c  (  (ojUjy  +  {opnjy 


(af  n,)')  ^ 


(ojTijy  +  (opnjY  +  [opnjY  =  1  (71) 


which  shows  that  the  square  of  the  acoustic  eigenvalues  (68)  equals  the  sum  of  the  square  of  the 
streamline  and  crossflow  acoustic  eigenvalues  (70).  On  any  plane  11  that  contains  a,  and  n, 
therefore,  acoustic  propagation  naturally  decomposes  into  two  directional  components,  one  along 
a  and  the  other  along  as  illustrated  in  in  Figure  4. 
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Figure  4:  Polar  Variation  of  Square  of  Acoustic  Speeds  on  Plane  11 


This  planar  two-direction  acoustic  decomposition  remains  unaltered  on  plane  11,  as  this  plane 
spans  the  entire  3-D  space  by  rotating  about  the  line  of  a.  Acoustic  propagation,  therefore,  decom¬ 
poses  into  a  directional  component  along  a  and  an  isotropic  component  on  any  plane  orthogonal 
to  a.  This  result  also  follows  from  the  expression 


(of  rij)  ^  -I-  (af  ny) ^  =  1  -  (ujUj)  ^  (72) 

from  (71),  which  shows  that  acoustic  propagation  on  the  plane  of  o^i  and  only  depends  upon 
(ajUj),  hence  on  the  angle  between  a  and  n.  For  any  such  angle,  the  “Ihs”  of  (72)  remains  constant 
for  any  orientation  of  n  with  respect  to  and  This  result  indicates  that  (72)  is  the  equation 
of  a  circle,  which  signifies  isotropic  acoustic  propagation  on  the  plane  of  and  . 

For  a,  ,  and  0^2  respectively  pointing  in  the  streamline  and  two  mutually  perpendicular 
crossfiow  directions,  the  Euler  fiux  divergence  can  then  be  decomposed  as 


^0^  =  aA^ajUk 


(73) 


The  weights  1  and  a,  1  <  a  <  0,  respectively  for  the  crossflow  and  streamline  components  in  this 
expression  are  difierent  from  each  other  because  the  streamline  and  crossflow  characteristic  velocity 
components  remain  distinct  from  each  other,  following  the  Euler  eigenvalues  (42);  furthermore, 
the  pressure-gradient  term  in  decomposition  (48)  will  contribute  difierent  characteristic  velocity 
components  along  the  streamline  and  crossflow  directions.  The  magnitudes  of  needed  acoustic 
upstream  bias  for  (73)  along  these  three  directions,  therefore,  will  have  to  differ  from  one  another 
and  this  differential  upstream  bias  can  be  easily  instituted  through  the  distinct  weights  a  and  1 
on  the  streamline  and  crossflow  components. 
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5.3.2  Absolute  Acoustics  Matrices 

Despite  its  zero  eigenvalues,  AjOj  features  a  complete  set  of  eigenvectors  X,  detailed  in  Appendix 
C,  and  thus  possesses  the  similarity-transformation  form 


A^aj  =  XA“X-^  =  A'A“+ -f  AA“-  ,  A“  =  A“+  -I-  A“- 


(74) 


where  the  diagonal  matrices  A“+  and  A“-  only  contain  ±c  along  their  diagonals.  The  matrices 
XA°‘+X~^  and  AA“-A“^  respectively  correspond  to  the  “forward”  and  “backward”  acoustic- 
propagation  matrix  components  of  AjCj.  A  bi-modal  upstream-bias  approximation  of  AjOj,  there¬ 
fore,  readily  follows  from  instituting  a  forward  and  a  backward  upstream-bias  approximation  re¬ 
spectively  for  the  forward-  and  backward-  propagation  matrices  in  (74).  Results  similar  to  (74) 
then  readily  follow  by  replacing  a  with  and  This  bi-modal  approximation  directly  yields 
both  the  non-negative-eigenvalue  matrices  U“aj  =  X  (A“+  —  A“") 


=  X„,(A“+-A‘‘-)X 


-1 


-1 


(75) 


and  the  associated  matrix  product  AjOj  a^dq/dx}-.  The  matrices  in  (75)  correspond  to  the  stream¬ 
line  and  crossflow  absolute  acoustics  matrices.  With  reference  to  the  results  in  Appendic  C,  the 


unabridged  form  of  the  matrix  product 


A“aj 


akdqfdxk  has  been  determined  as 


lAjOj  ak 


dxk 


pjd 
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,  a\c  , 

d^d^c  , 
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0 

0 

,  Q2dlC  , 

ajc  , 

d2(lsC  , 

0 

0 

,  azdic  , 

ajc  , 

0 

V  (c^-Pp)Pp/(cP£) 

,  0  , 

0  , 

0 

(c^-P.)  /c  ) 

dk 


dq 

dxk 


(76) 

with  analogous  results  obtained  by  replacing  a  with  and  0^2.  This  result  is  considerably 
simplified  upon  imposing  on  it  the  conditions  that  its  continuity  and  energy  components  satisfy 
the  acoustic-field  results  (59)  and  that  its  momentum  components  satisfy  results  (62)  and  (63)  on 
the  basis  of  a  wave-like  field  approximation  outlined  below. 

Upon  imposing  that  (59)  should  also  apply  to  (76)  for  vanishing  Mach  numbers,  hence  negligible 
flow  kinetic  energy,  the  continuity  and  energy  components  in  (76)  become 


Po  dp  Pe  be  _  dp 

dk  n  "I"  o  —  o 

c  dxk  c  dxk  oxk 


cPe 


P,  dp  ,  dE  dE 

dk—  +  - -  (^k-^  =  cak^  (77) 


dxk 


dxk 


dxk 


Next,  imposing  that  (62)  and  (63)  should  also  apply  to  (76)  under  the  same  conditions,  simplifies 
the  momentum  components  in  (76).  For  instance,  the  x\  momentum  component  in  (76)  becomes 


2  dmi  dm2  dmz 

CO-^Qk  o  CCl\0/2dk  “I"  CUidzdk' 


cafok 


dxk 
dmi 


dxk 


dxk 


2  dmk  2  dmk  dmi 

+  caiaa-^:—  +  caia^-^—  =  cafc- 


dxk  '  dxk  '  dxk  dxk 

Analogous  results  are  obtained  for  the  X2  and  xz  momentum  Components,  hence 

dmj  _ 


dmi 

cxiiCLjd}^-^  — 

dxk  oxk 


1  <  z  <  3 


(78) 


(79) 


Following  (62),  these  results  are  exact  for  wave-like  acoustic  fields.  For  non-acoustic  fields,  the 
wave-like  field  approximation  for  (78)  rests  upon  the  first  term  of  a  mean- value-theorem  expansion 
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of  a  general  solution  q  =  q{xi,X2,xs,t).  Appendix  A  shows  that  such  a  solution  can  be  expressed 
as  g  =  where  is  expressed  as  in  (38).  A  mean-value-theorem  expansion  of  this 

general  solution  thus  depends  upon  jys,  and  774  and  yields 

,  0<  62,93,64  <1  (80) 

where  772, 773,  and  774  are  now  fixed,  which  render  9(771,  93, 774)  a  wave-like  function  with  respect  to 

771,  for  which  (62)  and  (63)  apply.  The  wave- like  field  approximation  also  rests  on  the  recognition 
that  the  eigenvalues  of  the  absolute  acoustics  matrix  A“aj  achieve  their  maximum  along  the 
streamline  direction  a,  as  shown  in  Section  3,  which  justifies  the  choice  n  =  a  in  the  expression 
771  =  XjUj  —  Xt  for  expansion  (80). 

Relations  (77)-(78),  therefore,  lead  to  the  beautifully  simple  acoustic-field  result 


q{‘nu‘n2,ri3,m)  =  q{vi^m,fj3,rj4)  +  Y^^AT]e 


dq  _  T  ^9  _  dq 

OXk  OXk  OXk 


(81) 


A  simplification  similar  to  (77)-(80)  is  also  applied  to  the  crossflow  components  ^  and 

leads  to  the  similar  result 


Ni  dq 


Ni  dq 


For  the  stated  approximations,  therefore,  results  (81)-  (82)  indicate  within  matrix  product  (76) 
and  its  crossflow  counterpart  the  equivalence  of  replacing  A“aj|,  and  with  the 

matrix  cl,  of  which  all  eigenvalues  approach  +c  for  vanishing  Mach  number.  Appendix  C  also 
presents  a  different  choice  for  A“+  and  A“-  that  makes  the  absolute  acoustic  matrices  identical  to 
cl  without  any  need  for  resolutions  (77)-(80). 

For  acoustic  flows  and  related  dependent  variables  p,  mi,  m2,  m3  and  E  expressions  (81)- 
(82),  based  on  (74),  axe  exact.  Figure  5  in  Section  (7.4)  then  reveals  that  the  contribution  to  the 
upstream-bias  formulation  associated  with  (82)  remains  significant  for  M  <  1.0  only,  whereas  the 
contribution  associated  with  (81)  remains  significant  for  M  <  0.39  only.  These  computationally 
advantageous  approximation  on  the  acoustics  upstream-bias,  not  on  the  flux  divergence  itself, 
therefore,  will  be  used  essentially  in  these  Mach  number  ranges. 


6  Acoustics-Convection  Flux  Divergence  Decomposition 

The  acoustics-convection  flux  Jacobian  decomposition  consists  of  components  that  genuinely  model 
the  physics  of  multi-dimensional  acoustics  and  convection.  These  components  combine  the  com¬ 
putational  simplicity  of  FVS  with  the  accuracy  and  stability  of  FDS  and  also  feature  eigenvalues 
with  uniform  algebraic  sign.  This  formulation  eliminates  the  unstable  linear-dependence  problem  in 
steady  low-Mach-number  flows  and  satisfies  by  design  the  upstream-bias  stability  condition.  As  the 
Mach  number  increases,  the  formulation  smoothly  approaches  and  then  becomes  an  upstream-bias 
approximation  of  the  entire  flux  divergence,  along  one  single  direction. 
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6.1  Equivalent  Decompositions 

The  previous  sections  have  shown  that  the  Euler  flux  divergence  can  be  equivalently  expressed  as 


dxi 


3  _ 


dx 


dx  n 
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(1-/3) 
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dq 


=  < 
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+  K  Jf;,' + A“- Jf;,')  <=£+ 
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dq 


ag 


dxk 


dx 


dxk  ^  dxj 


(83) 

where  the  flrst  expression  is  convenient  for  a  characteristics-bias  approximation  within  the  stream¬ 
line  region,  for  high-subsonic  and  supersonic  Mach  numbers,  and  the  second  expression  is  conve¬ 
nient  for  a  characteristics-  bias  approximation  within  the  crossflow  region  and  within  the  streamline 
region,  for  low  subsonic  Mach  numbers. 

For  a  multidimensional  upstream-bias  approximation  throughout  the  flow  fleld  and  for  all  Mach 
numbers,  the  flux  divergence  can  thus  be  cast  as  the  linear  combination 


dx 
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J-  +  3^ 

dxj  ^  dxj 
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dfT 


I  -h  5  |a  (XA“+A-^  -H 


Ofc 


A. 
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(84) 


with  0  <  a,  5  <  1,  which  leads  to  the  following  acoustics-convection  flux  divergence  decomposition 

=  5a  (x A“+X-^  -f  XA“- 


dxk 


+5  (x„.  A‘*x;;  +  +  a  {x„^  A-*x;l  +  x„a‘-x;^)  a,"’  £-+ 


dq 


dxk 
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df^  df^ 


A  +  (1_5«)/3 


df^  df” 

— ' -I- (1 -5a)(l +  5(a- 1) — ~ 


dxj 


t-  5  (1  -  a) 


In  this  decomposition,  the  expression 


a/? 


'dxj"'^‘  '~'"^"^'~^dxk  dxj 

(85) 

is  counted  as  one  term  because,  with 


+  (1  -  5a)/34 

reference  to  (48),  the  eigenvalues  associated  with  this  term  will  all  keep  the  same  sign  within 
the  streamline  region,  since  (1  —  aa)l3  <  /3.  As  mentioned  in  Section  5.3,  an  upstream-bias 
approximation  for  the  Euler  flux  divergence  is  developed  by  establishing  upstream  approximations 
only  for  the  first  8  terms  in  (85).  One  justification  for  the  selective  upstream  formulation  on  these 
terms  rests  on  the  physical  acoustics  and  convection  significance  of  these  terms. 

Since  the  bi-modal  crossflow  wave  propagation  region  exists  for  all  Mach  numbers,  the  crossflow 
acoustic  term  a  A‘*X;;  +  A“-  X-‘)  4‘^+a  {x„^  A-*  Jf-'  +  X„^  A‘-  J?-' )  has 

to  receive  a  bi-modal  upstream  bias.  The  function  5  will  decrease  for  increasing  M,  Isecause  the 
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bi-modal  crossflow  region  shrinks  for  increasing  Mach  number.  For  supersoaic  flows,  the  flux  di- 

fa/?  a/?l 

vergence  ^  can  be  upstreamed  in  its  entirety  along  a  streamline,  which  implies  that  the 

weight  (1  —  aa)/3  must  approach  1,  as  M  increases  toward  M  >  1,  while  0  <  S,a,/3  <  1  simul¬ 
taneously.  On  the  other  hand,  the  weight  5  controls  the  magnitude  of  crossflow  upstream  bias, 
which  does  not  have  to  vanish  for  supersonic  flows,  because  of  the  presence  of  the  finite  bi-modal 
wave-propagation  regions.  Hence,  a  does  not  have  to  approach  0  for  these  flows,  which  implies 

that  (1  —  o;)  remains  less  than  1  for  finite  M  >  1.  As  a  consequence,  the  contribution  (1  —  3a)0^- 
is  the  only  correct  pressure-gradient  fraction  within  (85)  that  should  be  streamline  upstreamed  for 
subsonic  flows,  so  that  as  M  increases  the  pressure-gradient  weight  (1  -  5a) can  smoothly  and 
uniformly  approach  1,  with  0  <  5,a,/3  <  1. 

For  any  magnitude  of  both  pressure  and  pressure  gradient,  the  convection  field  uniformly  carries 
information  along  streamlines;  hence,  can  receive  an  upstream  bias  in  its  entirety  for  any 

M.  The  matrices  XA^+X-'^  and  XA'^-X-^  account  for  the  streamline  bi-modal  propagation 
of  acoustic  waves;  these  matrices  are  thus  used  for  an  acoustics  upstream-bias  approximation 
within  the  streamline  region,  for  low  Mach  numbers.  As  the  Mach  number  increases  from  zero, 
therefore,  a  smaller  and  smaller  fraction  5a  (XA“+X“^  -|-  ATA®- X“^)  of  (XA^+X"^  -f  XA“-X“^) 
is  upstream  approximated,  which  implies  that  5a  approaches  0  for  increasing  M,  and  5a  =  0  for 

>  1.  The  pressure  gradient  too  accounts  for  the  bi-modal  propagation  of  acoustic  waves, 

but  in  conjunction  with  As  the  Mach  number  increases  from  zero,  a  larger  and  larger  fraction 

(1  -  5a)/?^  of  the  pressure  gradient  can  thus  be  upstreamed  in  the  same  direction  as  and  along 
df^  df? 

while  (1  -5a)(l  -  P)-^  is  upstreamed  in  the  opposite  direction.  Hence,  the  function  /? 
has  to  increase  for  increasing  Mach  number.  It  follows,  therefore,  that  the  upstream-bias  functions 
5,  a  and  /3  have  to  depend  upon  M  and  ensure  physical  significance  of  the  overall  upstream-bias 
approximation  to  (85).  These  functions  are  then  cast  through  the  expressions 

S  =  {1  -  aa){2p  -  1)  ,  a  =  5a  ,  =  a  (86) 

These  substitutions  in  terms  of  the  functions  a,  a^  and  6  lead  to  significantly  simpler  expressions 
for  the  characteristics  fiux. 


6.2  Multidimensional  Characteristics  Euler  Flux 

With  reference  to  (19),  given  the  physical  significance  of  the  terms  in  decomposition  (85)  and 
algebraic  signs  of  the  corresponding  eigenvalues,  the  associated  principal  direction  unit  vectors  for 
these  terms  are 


ai  =  -02  =  07  =  -08  =  o  ,  03  =  -04  =  o'^i  ,  O5  =  -oe  =  0^2  ^  09  =  Oio  =  Oil  =  0 

(87) 

At  each  flow-field  point,  o  is  parallel  to  the  local  velocity  vector,  whereas  o^i  and  0^2  remain 
orthogonal  to  velocity. 

With  (85),  (87)  and  approximations  (81)-  (82),  the  general  upstream-bias  expression  (19)  di¬ 
rectly  yields  the  acoustics-convection  characteristics  flux  divergence 


dxi 


eip  (aoiOj  -t-  a"^  (af'aj"'  -h 


dXj  ^'dXj 


(88) 


In  this  result,  the  expressions  (caaiUj-^  Ui^  and  (ca"^ 

determine  the  upstream  biases  within  respectively  the  streamline  and  crossflow  wave  propagation 
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regions.  These  two  expressions  combined  then  induce  a  correct  upwind  bias  along  all  wave  propaga- 

a/p  . 

tion  regions.  In  particular,  the  coupling  of  an  upstream  approximation  for  (1  —  cc)/?^,  via  aj,  with 

d/? 

a  downstream  approximation  for  (1  -  a)(l  —  via  og  results  in  an  overall  upstream  approx¬ 

imation  of  the  pressure  gradient,  but  with  variable  weight  6.  The  operation  count  for  expression 
(88)  is  then  comparable  to  that  of  an  FVS  formulation.  The  terms  in  this  expression,  furthermore, 
directly  correspond  to  the  physics  of  acoustics  and  convection.  Expression  (88)  determines  ff 
itself,  up  to  an  arbitrary  divergence-free  vector,  as 


ff  =  Mq) 


etj) 


(^aaiUj  + 


NlNi 


+  af-a;" 


Jj  dXi 


+  d, 


df]  Off 

C/X  j  OfCj 


(89) 


According  to  this  result,  the  intrinsic  multi-dimensionality  of  each  component  ff  derives  from  its 
dependence  upon  the  entire  divergence  of  fj  and  /J. 

For  vanishing  Mach  numbers,  a  and  will  approach  1  whereas  5  will  approach  0.  Under  these 
conditions,  (88)  reduces  to 


dxi 


etj) 


'  dq  ^  a//' 

C-X - 1-  0,i— — 

I  uXi 


(90) 


which  essentially  induces  only  an  acoustics  upstream  bias.  Heed  that  this  bias  becomes  independent 
of  specific  propagation  directions,  for  it  no  longer  depends  on  ^aaiaj  +a^ 

This  bias,  therefore,  becomes  isotropic,  in  harmony  with  the  isotropic  propagation  of  acoustic 
waves.  Observe,  moreover,  that  the  components  within  dff  /dxj  remain  linearly  independent  of 
one  another,  which  avoids  the  linear-dependence  instability  in  the  steady  low-Mach-number  Euler 
equations.  For  supersonic  flows,  a  =  0  and  5  =  1  and  (88)  thus  becomes 


dxi 


etp 


+  d. 


^2 


dq 

dxj 


(91) 


which  depends  on  the  crossflow  component  of  the  absolute  acoustics  matrix  and  the  entire  diver¬ 
gence  of  the  Euler  inviscid  flux  vector. 


7  Infinite  Directional  Upstream  Bias 


In  Jacobian  form,  expression  (88)  becomes 


dxj  dxj 


dxi 


eip 


(aaiaj  +  -f  I  +  a., 


dq 


did 


dq  I  dxj 


(92) 


Expressions  (88),  (92)  essentially  depend  upon  the  six  upstream-bias  functions  01,02,03,0:,  5,  q^, 
considering  that  0^1  and  0^2  are  obtained  from  a  as  shown  in  Section  7.4.  In  order  to  ensure 
physical  significance  for  the  characteristics  flux  (89),  hence  for  the  upstream-bias  approximation 
to  decomposition  (85),  these  functions  are  determined  by  imposing  on  (92)  the  stringent  stability 
requirement  that  it  should  induce  an  upstream-bias  diffusion  not  just  along  the  principal  streamline 
and  crossflow  upstream  directions,  but  along  all  directions  n  =  (ni,  n2,  ns)  radiating  from  any  flow- 
field  point.  Appendix  D  shows  that  this  demanding  stability  conditions  is  satisfied  when  all  the 
eigenvalues  of  the  associated  upstream-bias  matrix 


A  =  ni\c 


-b  ^  d 


1  ^2  JV, 

+  «i  Oj 


)) 


dff  df^'' 

dq  dq  , 


Ui 


(93) 
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remain  positive  for  all  M  and  propagation  directions  n. 

Despite  the  formidable  non-linear  algebraic  complexity  of  A,  all  of  its  eigenvalues  have  been 
analytically  determined  exactly  in  closed  form.  Dividing  through  by  the  speed  of  sound  c,  the 
non-dimensional  form  of  these  eigenvalues  is 

■^0,2  =  ni  (auiOj  +  -t- rij  +  niaiVjUjM 

-^3,4  =  Tii  (auiCj  +  i^j+ 


+  S  (94) 

where  Vj  denotes  the  j***  direction  cosine  of  a  unit  vector  v  parallel  to  the  local  velocity  u.  Following 
the  considerations  after  (82),  in  particular,  all  these  eigenvalues  must  converge  to  1  for  vanishing 
Mach  number.  Heed  that  for  both  a  =  v  and  n  =  v,  the  functions  o;  and  5  within  (94)  determine 
the  corresponding  streamline  upstream-bias  eigenvalues 


Rather  than  prescribing  some  expressions  for  a  and  S  and  accepting  the  resulting  variations  for 
these  eigenvalues,  physically  consistent  expressions  for  the  streamline  upstream-bias  eigenvalues 
are  instead  prescribed  and  the  corresponding  functions  for  a  and  S  determined. 

7.1  Conditions  on  Upstream-Bias  Functions  and  Eigenvalues 

The  eigenvalues  (94)  are  expressed  as 


■^0,2  —  O'  +  M  ,  As, 4  =  a  -f 


(i+ V'’*) 


M±. 


,  1-5 


^  VjUjM  ± 


1-5 


p^VjUjM 


■^0,4  =  -^0,4  {M,n)  (96) 

to  stress  their  dependence  upon  both  M  and  n.  The  six  conditions  for  the  determination  of  the 
six  functions  oi,  as,  as,  o,  5,  are 

of  +  02  4-  03  =  1  ,  Ao,2  (M,  n)  >  0  ,  Ao,i  (M,  u)  =  Ai  ,  A4  (M,  w)  =  A4  ,  A3, 4  (M,  n)  >  0 

(97) 

where  Ai  and  A4  now  denote  prescribed  streamline  upstream-bias  eigenvalues.  The  first  condition 
stipulates  a  as  a  unit  vector  and  with  the  second  condition  it  determines  a,  and  for 
a,  ,  and  a^2  are  mutually  perpendicular  to  one  another.  In  particular,  these  two  conditions 
theoretically  confirm  that  the  unit  vector  a  has  to  point  in  the  streamline  direction,  which  implies 
that  a^i,  and  0^2  have  to  point  in  any  two  mutually  perpendicular  crossflow  directions.  The 
third,  fourth,  and  fifth  conditions  stipulate  that  the  streamhne  upstream-bias  eigenvalues  must 
equal  prescribed  eigenvalues,  and  thereby  determine  o;  and  5.  For  the  determined  a,  a^i,  a, 
and  5,  the  sixth  condition  then  establishes  a^. 

7.2  Streamline  Eigenvalue  A4 

This  eigenvalue  will  correlate  with  the  absolute  Euler  eigenvalue  \M  —  1|.  As  a  consequence,  A4 
will  vary  between  1  and  1  -  M  for  0  <  M  <  1  -  £m  and  smoothly  shift  from  1  -  M  to  M  - 1  within 
the  sonic  transition  layer  1  —  cm  <  M  <  1  -|-  where  £m  denotes  a  transition-layer  parameter; 


30 


Joe  lannelli 


in  this  work  em  =  5-  One  expression  for  A4  that  remains  smooth  and  meets  these  requirements  is 
the  composite  spline 


1-M  ,  0  <M<l-eM 

(M  —  1)^  Em  .  , 

A4(M)  =  <  — 2^ <  M  <  I  +  Em 

M-1  ,  1  +  em  <  M 


(98) 


7.3  Streamline  Eigenvalue  Ai 

This  eigenvalue  correlates  with  the  non-dimensional  Euler  eigenvalue  M,  but  it  too  has  to  equal  1 
for  M  =  0;  it  then  must  coincide  with  M  for  M  >  1  and  also  remain  greater  than  A4,  as  expressed 
through  (98),  for  consistency  with  the  Euler  eigenvalues  (42).  This  condition  in  particular  implies 
Ai  >  5.  It  thus  follows  that  Ai  will  vary  between  1  and  M  for  0  <  M  <  5  +sm-  An  expression  for 
Ai  =  Ai(M)  that  remains  smooth  and  meets  all  of  these  requirements  is  the  composite  spline 


(99) 


7.4  Upstream-Bias  Functions  a,  a^,  a,  6  and  cx^ 

These  functions  are  used  in  actual  computations  based  on  the  characteristics  flirx  divergence  (88). 
In  the  eigenvalues  Ao,2  in  (94),  the  components 

niaaiajTij  =  a  (ajUj)^  ,  niO.^ 

(100) 

are  already  non-negative  for  non-negative  o;  and  .  The  eigenvalues  Ao,2i  therefore,  will  remain 
non-negative  for  all  positive  a  and  a^,  including  a  0  and  >  0,  when  the  additional  com¬ 
ponent  niaiVjTijM  remains  non-negative  for  all  M.  This  requirement  is  met  along  with  the  first 
condition  in  (97)  when  a  =  v,  for 

niaiVjTijM  =  M  {vjnj)^  >  0  (101) 

This  finding  is  not  surprising,  for  the  streamline  direction  is  a  (  principal )  characteristic  direction. 
In  any  flow  region  where  velocity  vanishes,  it  is  then  computationally  convenient  to  set  ai  =  1, 
02  =  03  =  0.  The  components  of  a,  can  thus  be  expressed  as 


ui  U2  m 


(102) 
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It  is  straightforward  to  verify  that  with  these  expressions,  a,  ,  and  ^lre  unit  vectors  that 
remain  mutually  perpendicular  with  respect  to  one  another. 

Prom  Ai  and  A4  in  (95),  the  corresponding  expressions  for  both  a  =  a{M)  and  6  =  6{M)  are 
then  directly  and  exactly  determined  as 


^  (Ai(M)-A4(M))(Ai(M)-A4(M)+p„Af) 

a(M)  =  A.(M)-M  .  S(M)=^  (  A.(M)  -  A4(M)  ) 

where  according  to  the  fourth  and  fifth  conditions  in  (97)  the  streamline  eigenvalues  A4  and  Ai  are 
respectively  given  by  (98),  (99). 

Appendix  E  shows  that  the  sixth  condition  in  (97)  leads  to  the  following  expression  for 


a^{M)  = 


/ 3(cr^  (Mm)  -  1)  (Mm)  ^  (Mm)  (Mm)  -  1)  ^  ^^3 

Y  Mlf  Mm  )  \  Mlf  J 

0  <M  <  Mm‘, 

l(l  + - .  Mu<M 


where  superscript  prime  “ '  ”  denotes  differentiation  with  respect  to  M. 
The  variations  of  a  =  oi{M),  =  a^{M)  and  5  —  S{M)  in  Figure  5 


(104) 

indicate  that  these  three 


Figure  5:  Upstream-Bias  Functions 

functions  as  well  as  their  slopes  remain  continuous  for  all  Mach  numbers,  with  0  <  a, <  1 
and  O'  =  0  for  M  >  5  4-  em,  6  =  0  ior  M  <  0.3,  and  (5  =  1  for  M  >  1  -I-  em-  The  variation  of 
6  =  6{M)  shows  that  the  pressure-gradient  contribution  to  this  upstream-bias  formulation  increases 
monotonically,  while  remaining  less  than  25%  of  its  maximum,  for  0  <  M  <  0.7.  As  5  =  6{M) 
rises,  the  effect  of  the  acoustic-flow  result  (81)  is  concurrently  reduced,  because  the  streamline 
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upstream-bias  a  contribution  from  the  corresponding  acoustics  matrix  decreases  rapidly,  reducing 
by  75%  at  M  =  0.39. 

The  decrease  of  hence  of  the  crossflow  upstream-bias  from  (82),  is  less  rapid  because  this 
is  the  only  contribution  to  a  crossflow  upstream.  The  function  ,  nevertheless,  decreases  by  50%, 
at  the  sonic  state,  and  by  80%  for  M  =  1.8.  For  this  function,  forcing  non-negativity  of  A3, 4  as 
opposed  to  equality  to  a  prescribed  positive  constant,  in  particular,  ensures  minimal  cross  flow 
diffusion.  Expression  (104)  leads  to  the  conclusions 


lim  a^{M)  =  0 

M-^oo 


V  n 


(105) 


which  indicate  that  the  magnitude  of  crossflow  upstream  decreases  with  increasing  M.  This  result 
agrees  with  the  physics  of  high-M  flows,  where  the  bi-modal  propagation  region  narrows  about 
the  crossflow  direction,  as  seen  in  Section  4.4.  Convection  thereby  becomes  the  prevailing  wave 
propagation  mechanism,  which  therefore  reduces  the  need  for  acoustic  crossflow  upstream  bias. 


7.5  Streamline  Upstream-Bias  Eigenvalues 

As  Figure  6  shows,  the  streamline  upstream  bias  eigenvalues  (95)  remain  positive.  Furthermore, 


Mach  Number 


Figure  6:  Upstream-Bias  Eigenvalues 

these  eigenvalues  and  their  slopes  remain  continuous  for  all  Mach  numbers.  For  0  <  M  <  I  +  em^ 
the  eigenvalues  Ao,2)  A3,  A4  smoothly  approach  1  for  vanishing  M,  indicating  a  physically  consistent 
upstream-bias  approximation  of  the  acoustic  equations  (56)  embedded  within  the  Euler  equations. 
For  M  >  1  -|-  Em,  these  eigenvalues  respectively  coincide  with  the  Euler  flux  Jacobian  streamline 
eigenvalues  M,  M  +  1,  M  —  1,  which  corresponds  to  a  streamline  upstream-bias  approximation  of 
the  entire  flux  vector,  for  supersonic  flows.  Observe  also  the  smooth  transition  in  the  critical  sonic 
region  within  the  transition  layer,  where  A4  does  not  vanish,  but  remains  not  less  than  EmI‘^- 
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7.6  Polar  Variation  of  Upstream-Bias 

The  directional  variation  of  the  upstream  bias  eigenvalues  (94)  is  presented  in  Figures  7  -  8  for 
representative  subsonic  and  supersonic  Mach  numbers.  These  variations  are  obtained  for  a  variable 
unit  vector  n  =  (cos  0,  sin 0)  and  fixed  unit  vector  a  =  u,  in  this  representative  case  inclined  by 
+30°  with  respect  to  an  si  axis  on  any  plane  11  that  contains  the  velocity  vector. 


M  =  0.05  M  =  0.5 


Figure  7:  Polar  Variation  of  Subsonic  Upstream  Bias 


M=  1.5 


M  =  2.0 


Figure  8:  Polar  Variation  of  Supersonic  Upstream  Bias 
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These  figures  collectively  indicate  that  the  characteristics  flux  divergence  (88)  induces  a  phys¬ 
ically  consistent  upstream  bias  because,  for  any  Mach  number  and  wave  propagation  direction 
n,  the  associated  upstream-bias  eigenvalues  (94)  remain  positive  and  their  directional  variation 
mirrors  the  directional  variation  of  the  characteristic  Euler  eigenvalues  (42).  The  upstream-bias 
eigenvalues,  moreover,  are  symmetrical  about  the  crossflow  direction  and  characteristic  streamline, 
precisely  like  the  characteristic  Euler  eigenvalues  (10).  For  M  =  0.05,  the  directional  variation  of 
the  upstream-bias  eigenvalues  in  Figure  7  correlates  with  that  in  Figure  1  and  thereby  corresponds 
to  an  isotropic  upstream  bias,  in  complete  agreement  with  the  isotropic  acoustic  wave  propagation 
speed  in  the  Euler  equations.  For  increasing  Mach  numbers,  the  upstream  bias  becomes  anisotropic, 
again  in  agreement  with  the  anisotropic  distribution  of  the  Euler  eigenvalues  (42).  For  M  =  0.5 
this  anisotropy  is  already  evident  and  then  becomes  more  marked  for  supersonic  Mach  numbers 
as  indicated  in  Figure  8.  In  particular,  the  crossflow  upstream  bias  decreases  for  increasing  Mach 
number. 

Figures  9-10  compare  the  directional  variations  of  the  representative  upstream-bias  eigenvalue 
A3  and  the  corresponding  Euler  eigenvalue  A3  .  This  comparison  is  sufficient  to  depict  the  correla¬ 
tion  between  all  the  Euler  and  upstream-bias  eigenvalues,  for  Ao,2  and  Ao^2  topologically  similar 
to  each  other,  compare  Figures  2  and  8,  while  A4  and  A4  are  respectively  mirror  skew-symmetric 
to  A3  and  A3  with  respect  to  the  crossflow  direction. 


M  =  0.05 


110°  70° 


\ 


110° 


M  =  0.5 

100°  90°  S0° 


70° 


Figure  9:  Polar  Correlation  of  Subsonic  Characteristic  A3  and  Upstream  A3 
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Figure  10:  Polar  Correlation  of  Supersonic  Characteristic  A3  and  Upstream  A3 

As  Figures  9-10  indicate,  A3  is  symmetrical  about  the  characteristic  streamline,  precisely 
like  the  corresponding  characteristic  Euler  eigenvalue  A3  and  the  corresponding  polar  curve  is 
topologically  similar  to  the  Euler  eigenvalue  curve.  For  M  =  0.05,  A3  and  A3  virtually  coincide 
with  each  other  and  remain  direction  invariant,  which  corresponds  to  an  isotropic  upstream  bias 
in  correlation  with  the  acoustic  speed.  For  M  =  0.5,  Figure  9  indicates  that  A3  is  greater  than 
A3  in  the  streamline  direction.  This  results  corresponds  to  minimal  streamline  upstream  diffusion 
and  is  explained  in  terms  of  Figure  6,  which  shows  A3  in  the  streamline  direction  as  substantially 
less  than  M  +  1,  for  this  Mach  number. 

For  supersonic  Mach  numbers,  A3  in  the  streamline  direction  coincides  with  M  +  1.  As  shown  in 
Figme  10,  therefore,  the  magnitude  of  the  upstream  bias  for  supersonic  flows  is  virtually  identical 
to  the  magnitude  of  the  characteristic  eigen’values,  within  the  domain  of  dependence  and  range  of 
influence  of  any  flow  field  point.  Outside  this  region,  the  upstream-bias  eigenvalues  are  modestly 
less  than  the  characteristic  eigenvalues.  In  these  variations,  the  upstream-bias  eigenvalues  are 
vanishingly  small  in  the  cross-  flow  direction,  which,  in  particular,  corresponds  to  minimal  crossflow 
diffusion. 


8  Finite  Element  Weak  Statement 


The  divergence  (88)  of  the  characteristics  flux  leads  to  the  following  characteristics-bias  integral 
statement  with  implied  summation  on  repeated  subscript  indices  i,j 


Ju  [  dt  dxj  dxj 
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(106) 
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An  integration  by  parts  of  the  characteristics-bicis  expression  then  generates  the  weak  statement 


L 


W 


dq  dfj  dfv 

dt  dxj  dxj 


dfi 


.9f],  _  M 


+  Oj^ 

dxj 


_ l_ 

dxj 


dfi— 


c  (aa  •  naj  +  a'^  (a'^^  •  +  a  •  naj-^  +  a  •  najS-^^  j  dT  =  0 

(107) 

where  n  denotes  the  outward  pointing  unit  vector  orthogonal  to  the  boundary  dQ  of  Cl. 


8.1  Elimination  of  Spurious  Upstream  Boundary  Dissipation 

On  any  boundary  facet  where  w  ^0,  expression  (107)  would  induce  spurious  upstream  boundary 
dissipation  if  the  boundary  integral  were  not  evaluated.  It  is  possible,  however,  to  eliminate  both 
this  surface  integral  and  any  spurious  upstream  boundary  dissipation.  This  result  is  achieved  by 
setting  a,  a^i,  and  a^2,  within  every  finite  element  with  a  face  on  the  boundary  dCl,  equal  to  unit 
vectors  that  are  orthogonal  to  n.  One  choice  is 

{a)QQ  =  a  —  a-nn  ,  (a^i)^^  =  •  nn  ,  ■  nn  (108) 

With  this  choice  (o)an  •  n  =  0,  (0^1)90  •  n  =  0,  {a’^‘^)dn  •  n  =  0,  the  corresponding  surface  integral 
in  (107)  naturally  vanishes,  hence  no  spurious  upstream  boundary  dissipation  is  induced  by  not 
evaluating  this  integral. 


8.2  Galerkin  Finite  Element  Equations 

The  finite  element  weak  statement  associated  with  (107)  is 


(109) 


where  superscript  “h"  signifies  spatial  discrete  approximation.  The  approximation  exists  on  a 
partition  Cl^  C  Jl,  of  Cl.  This  partition  Cl^  has  its  boundary  nodes  on  the  boundary  dCl  of  Cl 
and  results  from  the  union  of  Ne  non-overlapping  elements  Cle,  Cl^  =  U^i  ^e-  Within  Cl^,  there 
exists  a  cluster  of  “master”  elements  Cl^ ,  each  comprising  only  those  adjacent  elements  that  share 
a  mesh  node  with  1  <  k  <  N,  where  N  denotes  the  total  number  of  mesh  nodes  and  hence 
master  elements. 

In  the  representative  case  of  a  2-D  flow,  as  Figure  11  shows,  the  discrete  test  function  within 
each  master  element  Cl^  will  coincide  with  the  “pyramid”  basis  function  Wk  =  Wk{x),  1  <  k  <  N, 
with  compact  support  on  Clj^.  Such  a  function  equals  one  at  node  Xk,  zero  at  all  other  mesh  nodes 
and  also  identically  vanishes  both  on  the  boundary  segments  of  Cl^  not  containing  Xk  and  on  the 
computational  domain  outside  Cl^. 
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Figure  11:  a)  Pyramid  Test  Function  for  fif ,  b)  Assembly  of  4  Uniform  Quadrilaterals 

The  discrete  solution  and  flux  at  each  time  t  assume  the  form  of  the  following  group 
linear  combinations 

N  N 

q''{x,t)  =  '^we{x)-c^{xi,t)  ,  fj{x,t)  =  '^we{x)-fj((^{xe,t))  (110) 

^=1  e=i 

of  time-dependent  nodal  solution  values  q^  {xe,  t),  to  be  determined,  and  trial  functions,  which  coin¬ 
cide  with  the  test  functions  wi{x)  for  a  Galerkin  formulation.  Similarly,  the  fluxes  /J  =  fj(q{x,  t)) 
and  /J  =  fj{q{x,t))  are  discretized  through  the  group  expressions 

h  ^  AT 

fj  =  ,  ff{x,t)  =  '^we{x)-ff(q'^{xi,t))  (111) 

e=i 

with  similar  expansions  applying  for  the  components  of  the  viscous  flux  /J^.  The  notation  for 
the  discrete  nodal  variable  and  fluxes  is  then  simplified  as  qe(t)  =  q(xe,t),  =  fj  (^q^{xi,t)^, 

^  expansions  (llO)-(lll)  are  then  inserted  into  (109), 

which  yields  the  discrete  finite  element  weak  statement 


dqe  ,  dw£ 
dt  dxi 


dfvh\ 
dXj  ) 


dfl  + 


dwkdwi  ft 


/n^  dxi  dx. 


eV*  [o*  (a‘a?a}  +  a"*  g,  +  aifl  +  ]  «  =  0 


for  1  <  k  <  AT.  There  are  three  implied  summations  with  respect  to  the  subscript  indices 
The  subscript  indices  i,j  in  this  expression  denote  cartesian-axis  directions,  hence  I  <  i,j  <  3, 
whereas  subscript  £  indicates  a  mesh  node,  hence  1  <  ^  <  iV. 

While  an  expansion  like  the  ones  in  (110)  for  V’'*,  a*,  (^,  a^,  a^i ,  a^2,  and  6^  can  be  directly 
accommodated  within  (112),  each  of  these  variables  in  this  study  has  been  set  equal  to  a  piece  wise 
constant  for  computational  simplicity,  one  centroidal  constant  value  per  element.  Likewise,  is 
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set  equal  to  a  reference  length  within  each  element,  typically  a  measure  of  the  element  size.  In 
this  study,  =  (^)e/2  within  each  element  “e”,  Where  {£)e  denotes  the  length  of  the  streamline 
diameter  of  the  generalized  ellipsoid  inscribed  within  the  element.  The  choice  of  a  diameter  in  the 
streamline  direction  for  (£)e  rests  on  the  recognition  that  the  streamline  is  a  characteristic  principal 
direction,  as  discussed  in  Sections  4.2  -  4.4. 

Since  the  test  and  trial  functions  we  are  prescribed  functions  of  x,  the  spatial  integrations  in 
(112)  are  directly  carried  out.  For  arbitrarily  shaped  elements,  these  integrations  take  place  via  the 
usual  finite  element  local-coordinate  transformation^^’^^  that  for  example  maps  a  quadrilateral 
into  a  square.  In  this  study,  the  resulting  coordinate-  transformation  metrics  within  each  element 
are  set  to  constants  equal  to  their  respective  centroidal  values.  This  simplification  allows  the  exact 
integration  of  the  remaining  integrals,  which  are  then  evaluated  only  once  for  each  computation. 
The  complete  integration  with  respect  to  x  transforms  (112)  into  a  system  of  continuum-time 
ordinary  differential  equations  (ODE)  for  determining  at  each  time  level  t  the  unknown  nodal 
values  1  <  £  <  N.  This  ODE  system  is  numerically  integrated  in  time  via  an  implicit 

diagonal  Runke-Kutta  algorithm  that  remains  absolutely  stable  for  stiff  non-linear  dissipative  sys¬ 
tems.  Concerning  the  boundary  variables,  no  extrapolation  of  variables  is  needed  in  this  algorithm 
on  a  variable  that  is  not  constrained  via  a  Dirichlet  boundary  condition.  In  this  case,  instead, 
the  finite  element  algorithm  (109)  naturally  generates  for  each  unconstrained  boundary  variable  a 
boundary-node  ordinary  differential  equation. 


8.3  Discrete  Conservation 


The  finite  element  equations  (112)  naturally  lead  to  a  discretely  conservative  algorithm,  for  they 
can  be  rearranged  in  a  traditional  “numerical  flux”  form.  When  the  mesh  nodes  of  a  rectangular 
grid  of  squares  are  numbered  via  integer  coordinates  i,  j,  for  a  representative  2-D  Euler  formulation, 
the  specific  numerical-flux  form  of  the  finite  element  system  is 


dq. 


dt 


+ 


F2.  1  —  F2.  ,  1 

Ax2., 


(113) 


which  is  analogous  to  the  form  of  upwind  finite- volume  schemes.^®  In  this  expression,  q.  .  denotes 
the  average  value  of  the  discrete  q^  at  node  The  finite  element  formulation  directly  generates 

this  average,  which  for  the  simple  representative  case  of  a  uniform  cartesian  mesh  with  constant 
Axi  and  Ax2  becomes 


+  Qi-i,j+i  +  ^Qi,3+i  +  g;+i.j+i 
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(114) 


Likewise,  the  formulation  directly  generates  the  intrinsically  multi-dimensional  “numerical  fluxes” 
jPi  ,  and  F2  ,  .  For  a  uniform  cartesian  mesh,  these  fluxes  are  expressed  as 


i+^,j 


Axe 


(115) 


and 
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—  ^e^caa2a£^^^  —  of  (116) 

».i+^  ij+i 

where  “(  )”  denotes  an  average  over  the  elements  sharing  node  Appendix  F  presents 

the  expressions  for  these  averages.  Apart  from  becoming  the  discrete  counterparts  of  the  multi¬ 
dimensional  characteristics-bias  flux  (19),  expressions  (115)-(116)  correspond  to  authentic  numer¬ 
ical  fluxes^®  because  they  respectively  reduce  to  fi{q)  and  f2{q)  when  each  nodal  value  of 
therein  is  replaced  by  q,  for  flxed  Axi  and  Ax2- 

The  non-negative  controller  ip  determines  the  amount  of  induced  upstream  bias  in  these  ex¬ 
pressions.  For  ^  =  0,  no  upstream  bias,  hence  dissipation,  is  induced  and  (109)-(116)  generate  a 
centered  approximation  for  the  Euler  equations.  The  local  element  length  hence  (£)e,  is  then 
chosen  so  that  ip  =  1  then  corresponds  to  a  fully  upwind  approximation,  hence  maximum  upstream 
dissipation.  The  numerical  solution  for  drives  the  local  magnitude  of  ip  within  each  element^^  : 
in  flow  regions  where  discontinuities  would  form  then  ip  =  V’max)  for  essentially  non-oscillatory 
shock  capturing,  and  in  smooth  flow  regions  then  ip  =  i/’min?  for  increased  accuracy.  Equations 
(109)-(116)  induce  an  upstream  bias  in  the  correct  direction  with  respect  to  the  velocity  direction 
because  the  direction  cosines  oi  and  02  of  velocity  u  determine  the  upstream  direction.  In  par¬ 
ticular,  when  either  oi  or  02  vanishes,  then  either  Fi  or  F2  becomes  a  centered  numerical  flux  in 
the  xi  or  X2  cartesian  direction,  while  featuring  an  acoustics  upstream  bias  in  the  corresponding 
orthogonal  direction,  as  induced  by  either  02  or  of.  In  any  case,  equations  (109)-(116)  ultimately 
induce  an  upstream  bias  on  the  entire  Euler  flux  divergence  for  any  orientation  of  a  =  w. 


9  Discrete  Upstream-Bias  Controller  -0^ 

This  section  introduces  a  new  upstream-bias  controller  ip^.  This  controller  varies  in  the  range 
0  <  ipiain  <  ^  V’max  <  1  and  controls  within  (112)  the  amount  of  induced  upstream-bias,  hence 

artificial  diffusion.  By  design,  ip^  =  0  corresponds  to  a  classical  centered  discretization,  hence  no 
upwind  diffusion;  ip^  =:  1,  instead,  corresponds  to  a  fully  upwind  formulation,  hence  maximum 
diffusion. 

Denote  then  with  ipi  the  numerical  value  of  the  controller  at  the  representative  node  “i” ,  for 
a  1-D  formulation,  or  at  the  face  shared  by  two  2-D  or  3-D  elements.  By  analogy  with  (110),  the 
discrete  controller  ip^{s,t)  is  cast  eis  the  following  linear  expansion 


2  N 

=  ^u;j(s)V’(sj,t)  =  Y^Wj{s)ipj  (117) 

i=i  j=i 


This  is  a  1-D  expansion  even  for  a  3-D  flow  and  corresponds  to  the  variation  of  ip^{s,t)  within 
each  element  with  respect  to  the  coordinate  s  along  the  arc  of  streamline  that  crosses  the  element 
centroid  and  intersects  two  faces  of  the  element;  the  centroidal  evaluation  of  this  expression  within 
each  element  then  yields 


,  +  tpi+i 

M  =  — 2 — 


(118) 


In  regions  of  smooth  flow  ip^  approaches  ipmin  for  a  local  reduction  of  upstream-bias  diffusion;  in 
region  of  discontinuous  solution  ip^  approaches  i/’max?  for  an  essentially  non-oscillatory  resolution  of 
local  discontinuities.  The  controller  will  thus  correlate  with  a  local  measure  (pi  of  slope  discontinuity. 
For  a  smooth  variation,  ipi  is  expressed  as  a  function  of  the  measure  (pi  as  the  composite  spline 
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V^max  "t" 


V^min 

V^max  “  '0niin 


=  { 


i^D 

-  {<Pd  -  Pi  +  ^  (fD  +  Pc)  Pi  -  2<^i 

V’max 


)  Pi  <  Pc 

,  Pc  <  Pi  <  Pd 
^  Pd  ^  Pi 


(119) 

where  0  <  v^j  <  1  and  the  positive  and  (pj^  respectively  denote  threshold  continuity  and 
discontinuity  measures.  Figure  12  shows  the  smooth  variation  of  this  type  of  spline  controller. 


Figure  12:  Variation  of  Controller  tpi 

The  implementation  of  ip^  thus  requires  a  set  of  points  within  where  the  slopes  of  the 
approximate  solution  are  generally  discontinuous.  For  the  finite  element  approximation  (110),  this 
set  of  points  is  taken  as  the  the  set  of  locations  on  faces  sheired  by  any  two  distinct  elements  in 
for  the  continuous  expansion  (110)  changes  approximating  polynomial  from  element  to  element, 
which  implies  that  solution  slopes  are  generally  discontinuous  at  element  faces.  This  type  of  slope 
discontinuity  is  depicted  in  Figure  13  via  the  local  normal  unit  vectors  and  respectively 
to  the  left  and  right  of  the  slope-discontinuity  point  P,  where  there  exist  two  distinct  normal 
unit  vectors,  one  for  each  of  two  elements  sharing  P.  This  picture  refers  to  1-D,  2-D  and  3-D 
formulations.  In  1-D,  point  P  in  this  picture  corresponds  to  a  finite  element  node;  in  2-D  and  3-D, 
point  P  is  the  centroid  on  the  face  shared  by  any  two  contiguous  elements  and  the  corresponding 
arcs  indicate  the  variation  of  along  the  element  centroidal  streamline  arc. 

With  reference  to  Figure  13,  the  magnitude  ||n^  — n^||  of  the  vector  difference  n^—n^  becomes 
proportional  to  a  bounded  measure  of  local  slope  discontinuity.  If  the  graph  slope  is  continuous 
at  P,  then  coincides  with  and  ||n^  -  ti^||  vanishes.  On  the  other  hand,  when  a  slope 
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X 


Figure  13:  Slope  Discontinuities  and  Local  Unit  Vectors 


discontinuity  exists  at  P,  as  shown  in  the  figure,  ||n^  —  n^||  varies  between  0  and  2  depending 
on  the  magnitude  of  the  slope  jump.  A  positive  measure  of  slope  discontinuity  that  vanishes  for 
continuous  slopes,  remains  bounded  and  strictly  varies  between  0  and  1  can  thus  be  defined  as 

<Pi  =  2 11*^  ~  ^  IU=a:i  (120) 

By  virtue  of  the  law  of  cosines,  the  local  measure  ipi  also  equals 

(121) 

where  6  denotes  the  angle  between  the  unit  vectors  in  Figme  13. 

With  reference  to  Figure  14  and  (121),  specific  numerical  values  for  ipj^,  ^min  aJid  V’max 
can  be  easily  established.  At  a  point  of  solution  smoothness,  like  point  i  —  1  in  the  figure, 
will  be  parallel  to  n^,  hence  6  =  0°  which  from  (121)  leads  to  =0.  At  a  shock,  instead,  9 
can  become  greater  than  90°,  as  shown  in  the  figure  for  point  i.  The  threshold  9  =  90°  is  thus 
selected  for  which  from  (121)  leads  to  =  l/\/2.  Numerous  numerical  experiments  with  the 
acoustics-convection  algorithm  have  indicated  that  a  minimal  amount  of  “backgroimd”  upstream 
bias  is  necessary  for  convergence;  this  finding  is  not  surprising,  since  the  formulation  is  essentially 
centered,  hence  devoid  of  any  upstream-bias  dissipation  for  =  0.  Hence,  i/’min  >  0,  with  typical 
numerical  values  in  the  range  |  <  V’min  <  5-  Concerning  V’max)  a  relation  with  V’min  readily  follows 
from  the  requirement  that  in  the  neighborhood  of  a  shock  the  maximum  upstream  bias  can  at 
most  correspond  to  a  fully  upwind  algorithm,  for  an  essentially  non-oscillatory  capturing  of  shock 
waves.  Hence,  from  (118),  V’i+i  <  1  with  i  +  \  denoting  the  centroid  of  a  finite  element  (  cell  ) 
that  supports  a  shock.  For  a  typical  case  of  a  shock  captured  within  at  least  two  cells,  as  shown 
in  Figure  14,  (119)  leads  to  ipi  =  V’max  and  ~  V’min-  From  (118),  therefore 


'0max  "b 


V’max  ^  2  '01) 


which  linearly  decreases  as  a  function  of  increasing  0inin-  The  specific  objective  of  letting  0^  vary 
as  the  solution  evolves  is  to  minimize  induced  upstream-bias  dissipation  for  maximum  accuracy 
within  the  prescribed  computational  stencil.  As  its  distinguishing  design  feature,  the  acoustics- 
convection  upstream  resolution  algorithm  nevertheless  remains  an  authentic  characteristics-bias 
formulation  for  any  0^  with  0inin  <  0^  <  0max- 
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Figure  14:  Local  Unit  Vectors  at  a  Shock 


The  general  expression  of  ipi  corresponding  to  a  scalar  component  of  directly  derives 
from  the  finite  element  expansion  (110),  which  can  be  expressed  in  synthetic  implicit  form  as 
F{q^,x,t)  =  q^  —  q^{x,t)  =  0.  Hence,  a  normal  unit  vector  n  can  be  cast  at  each  time  level  t 
as  n  =  gradF'(9g,®,<)/||gradF||,  where  the  vector  operator  “grad”  encompasses  the  dependent 
variable  q^.  For  an  n-dimensional  formulation,  with  1  <  n  <  3,  the  expression  for  the  corresponding 
(fi  at  time  level  t  and  face  centroid  Xi  then  becomes 


ipi  =  {Xi,  t) 


The  1-D  form  of  (pi  is 


1 


dxt  dxt  \\ 


1  ^  dqc^  dq^^ 

dxi  dxi  j 


+  E 

i=i 


dqc^ 
dx  <1 


2  n 


dq: 


hL 


dxj 


(123) 


X=Xi 


1 


[  dx  )  J 
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dq^c^ 

dx 

dq'^^  ^ 

dx 

2  1 

(124) 


where  the  partial  derivatives  are  determined  through  the  finite  element  expansion  (110),  and  where 
superscripts  L  and  R  indicate  evaluation  within  the  elements  respectively  to  the  left  and  right  of 
the  shared  face  where  x  =  Xi  lies. 

The  representative  1-D  form  of  (124)  at  node  i  of  a  uniform  grid  is 


<Pi 


+ 


Ax 


Ax 


[V  V^Aa:2  +  (g 

Ct+1  Qcif  \l +  {qci  -  qci.iY 
gcj+l  ~  gcj _  _ gcj  ~  gcj-l _ 


y'Ax2  +  (9ci+i  -fc)^  \/Aa:2  +  (gc,  -gci_i)' 


(125) 


where  the  denominator  never  vanishes.  This  expression,  furthermore,  remains  bounded  and  dif¬ 
ferentiable  for  arbitrary  nodal  values  of  q^.  For  the  sole  purpose  of  determining  the  order  of  this 
expression  with  respect  to  Ax,  for  a  smooth  solution  over  any  two  contiguous  elements,  the  discrete 
solution  values  ^Cj ,  ^  ~  1  <  i  <  *  + 1)  over  these  elements,  can  be  considered  as  the  nodal  values  of  a 
single  auxiliary  continuous  functions  qc{x,t):  a  Lagrangian,  trigonometric,  or  other  interpolant  of 
the  qc/s  over  both  elements.  With  this  consideration,  the  Taylor’s  series  expansion  of  (125)  yields 


V’i  = 


9c 


ft 

c 


Ax 


1+  (9c(®i><))' 


+  O  (Ax^)  =  /C  1  -h  {q'cixi,t)y  ^  ^  +  ^ 


(126) 


where  superscript  prime  indicates  differentiation  with  respect  to  x  and  K  denotes  the  local  curva¬ 
ture.  This  expansion  reveals  that  tpi  decreases  for  vanishing  Ax.  Even  for  large  slopes,  furthermore, 
ipi  remains  of  order  Ax  in  regions  of  small  curvature.  Only  when  both  curvature  and  slope  drasti¬ 
cally  rise,  e.g.  at  a  shock,  will  ipi  increase,  which  precisely  corresponds  to  the  desired  behavior. 


10  Implicit  Runge  Kutta  Time  Integration 

The  finite  element  equations  (112)  along  with  appropriate  boundary  equations  and  conditions  can 
be  abridged  as  the  non-linear  ODE  system 

M^^  =  F{t,Q{t))  (127) 

where  corresponds  to  the  coupling  of  time  derivatives  in  (112),  and  F  {t,Q{t))  represents 

the  remaining  terms  in  (112).  The  numerical  time  integration  of  (127)  in  this  study  takes  place 
through  a  new  class  of  two-stage  diagonally  implicit  Runge-Kutta  algorithms^  (IRK2)  expressed 
as 


Qn+l  Qn 
MKi 

MK2 


biKi  -b  62-^2 

At  •  F  {tfi  +  Cl  At,  Qn  +  oii.R^i) 

At  •  F  {tn  +  C2At,  Qn  +  a2i-R^i  +  022^2) 


(128) 
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where  n  now  denotes  a  discrete  time  station  and  bi,  62,  ci,  C2,  an,  021,  and  022  indicate  constant 
Runge-Kutta  coefficients,  subject  to  the  constraints  c,  =  bi  =  1.  The  coefficients 

for  second  order  accuracy  are  listed  in  the  Table  1.  With  these  coefficients,  in  particular,  algorithm 

Table  1:  Runge-Kutta  Coefficients 


61  62  an  a2i  022 

IRK2 

3-x/3  l-^-^/3  3-V3  ^  v/3-1 

4  4  6  ^  2 

(128)  becomes  absolutely  stable  for  arbitrary  stiflF  non-linear  dissipative  ODE  systemsJ’^^ 

This  algorithm  is  implicit  because  the  entries  in  the  arrays  Ki  and  K2  remain  coupled  and  are 
then  computed  by  solving  algebraic  systems.  Diagonally  implicit  signifies  that  Ki  is  determined 
independently  of  K2-  Thus,  given  the  solution  Qn  at  time  Ki  is  computed  first,  followed  by 
K2.  The  solution  Qn+i  is  then  determined  by  way  of  the  first  expression  in  (128). 

The  terminal  numerical  solution  is  then  determined  using  Newton’s  method,  which  for  the 
implicit  fully-coupled  computation  of  the  IRK2  arrays  Ki,  1  <  i  <  2,  is  cast  as 


M.  —  an  At 


Kf)  =  At  F  {tn  +  Ci  At,  Of)  -  MKf 


QP  =  Qn  +  anKP  +  ai2KP  (129) 

where  =  0  for  j  >  i,  p  is  the  iteration  index,  and  Kf  =  Ki  for  i  =  2;  for  linear  finite  elements, 
the  Jacobian 

/dF\‘P 

Ji(Q)  =  M-anAt(^—J^^  (130) 

then  becomes  a  block  sparse  matrix.  The  initial  estimate  Kf  can  be  set  equal  to  the  zero  array. 
When  one  iteration  only  is  executed  for  (129)  within  each  time  interval  Newton’s  iteration  becomes 
akin  to  a  classical  direct  linearized  implicit  solver. 


11  Concluding  Remarks 

The  characteristics-bias  upstream  method  rests  upon  the  physics  and  mathematics  of  multi  -  di¬ 
mensional  characteristic  acoustics  and  convection  for  general  equations  of  state.  It  generates  the 
upstream  bias  at  the  difierential  equation  level,  before  any  discrete  approximation,  by  way  of  a 
characteristics-bias  system  and  associated  decomposition  of  the  Euler  fiux  divergence  into  convec¬ 
tion  and  streamline  and  crossflow  acoustic  components. 

A  natural  Galerkin  finite  element  discretization  of  the  characteristics-bias  system  directly  pro¬ 
vides  an  accurate  and  genuinely  multi-dimensional  upstream-  bias  approximation  of  the  Euler  and 
Navier-Stokes  equations,  in  the  form  of  a  non-linear  combination  of  upstream  diffusive  and  down¬ 
stream  anti-diffusive  flux  differences,  with  greater  bias  on  the  upstream  diffusive  flux  difference. 

Along  all  the  infinite  directions  of  wave  propagation,  the  formulation  induces  anisotropic  and 
variable-strength  consistent  upstream  bias  that  correlates  with  the  spatial  distribution  of  charac¬ 
teristic  velocities.  The  developments  in  this  report  have  implemented  the  algorithm  using  a  linear 
approximation  of  fluxes  within  quadrilateral  cells  without  any  MUSCL-type  local  extrapolation  of 
variables. 
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This  characteristics-bias  algorithm,  admits  a  straightforward  implicit  implementation,  features 
a  computational  simplicity  that  parallels  a  traditional  centered  discretization,  and  rationally  elim¬ 
inates  superfluous  artificial  difliision. 
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A  Solutions  in  Non-linear  Wave-Like  Form 

With  implied  summation  on  repeated  indices,  a  solution  of  the  Euler  equations  in  non-linear  wave¬ 
like  form  is  expressed  as 


9  =  Q(m)  ,  m=x-n-  X{q)t  =  XjUj  -  X{q)t  (131) 

where  n  denotes  a  propagation-direction  unit  vector,  independent  of  {x,t),  and  A  =  X{q)  indicates 
a  solution-dependent  wave-propagation  velocity  component  along  the  n  direction.  For  x-n  —  X{q)t 
equal  to  a  constant,  the  wave-like  solution  q  in  (131)  remains  itself  constant,  which  also  makes 
X  =  X{q)  a  constant.  In  the  time-spax;e  continuum  {x,t),  therefore,  the  equation 


x-n-Xt  =  C  (132) 

represents  a  hyper-plane  with C  and  (ni, n2, n3,—X(q))  respectively  denoting  a  constant  and  a  unit 
vector  orthogonal  to  this  plane,  as  depicted  in  Figure  15.  The  surface  that  remains  tangent  to 
these  planes  for  all  n ’s  is  then  tangent  itself  to  a  characteristic  surface,  for  A  corresponding  to  a 
characteristic  velocity  component.  Equation  (132)  corresponds  in  the  (xi,X2)  space  of  a  2-D  flow 
to  a  bundle  of  parallel  lines  for  a  fixed  t  and  dilferent  C.  For  a  fixed  C  and  variable  t,  the  equation 
represents  a  C-line  propagating  in  the  n  direction  with  velocity  component  A. 

Next,  consider  the  coordinate  transformation 


7?1  =  XjUj  -  x{q)t  ,  7/2  =  »?2  iXl,X2,  X3,t)  ,  7/3  =  7/3  (xi,  X2,  X3,  t)  ,  7/4  =  7/4  (Xl,  X2,  X3,  t) 

(133) 

with  7/2,  7/3,  and  7/4  chosen  so  that  the  reference  frame  (7/1,7/2,7/3,774)  is  also  right-handed  like 
(xi,X2,X4,t).  With  this  transformation,  q{xjnj  -  X{q)t)  =  9(7/1)  and  for  each  t  and  non-linear  g, 
the  7/1  axis  then  points  in  the  (711,712,713, —A(9))  direction.  This  result  follows  from  the  partial 
derivatives  of  7/1  with  respect  to  (a;i,X2,X3,t) 


-Hq)  , 

^1  +  t 

9A 

dn\. 

\  9m 

)  dXj 

(134) 

-A 

Til 

n2 

n3 

(135) 

1 

1 

dm  ~~ 
dxi 

dm 

dX2 

dm 

dX3 

which  shows  that  the  7/1-axis  direction  vector  (|^,  ^)  is  parallel  to  (711,712,713,  -X{q)). 

Since  ti  is  a  unit  vector,  the  second  result  in  (134)  also  leads  to  the  following  relations 


97/1 

axi 


1 


97/1  97/1 

dxj  dx£^^^^ 


(136) 
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Figure  15:  Wave  Plane  and  Reference  Frames 


The  partial  derivatives  of  the  wave-like  solutions  (131)  also  satisfy  the  conditions 


w  X  9g 


w  X  dq  dq 


This  result  follows  from  the  partial  derivatives  of  (131)  as 


dt 


dg  _  dq  f  dX  dq 
dxj  dqi  y  dq  dxj 


(137) 


(138) 


where  ^  denotes  a  column  array  and  ^  indicates  a  row  array,  which  implies  that  is  a 

square  matrix.  These  relations  then  lead  to  the  two  systems 


dq  5A\ 
dm  dq) 


dt 


dq  9A\ 
dm  dq) 


dq  _  dq 
dxj  dm 


(139) 


where  I  denotes  a  4  x  4  identity  matrix.  Since  the  square  matrix  ^  ^  results  from  the  product 
of  column  and  row  arrays,  the  rows  of  this  matrix  are  all  linear  combinations  of  the  single  ^  row. 
All  the  eigenvalues  of  this  matrix,  but  one,  therefore,  vanish.  Consequently,  all  the  eigenvalues  of 
the  matrix  (l  +  equal  1,  except  one  eigenvalue  A,  which  has  been  exactly  determined  as 


dm 


(140) 


As  the  product  of  all  its  matrix  eigenvalues,  therefore,  the  determinant  of  the  matrix  ^7  + 

equals  A,  and  the  associated  matrix  of  cofactors  is  then  denotes  with  B.  Systems  (139)  can  then 
be  solved  to  yield 


A(g)  „  dq  dq  ^  nj  dq 
dm  ’  dxj  l+<^  dm 


dt 


(141) 
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With  a  finite  if  the  denominator  1+^^  approaches  zero,  then  ||  and  ^  become  unbounded, 

which  indicates  that  a  positive  t  —  ~1/^  corresponds  to  a  “breaking”  time  of  shock- wave  forma¬ 
tion.  In  any  region  that  does  not  contain  shock  waves,  or  in  the  absence  of  shock  waves,  therefore, 
A  =  1  4- 1-^  does  not  vanish,  the  matrix  (j  +  is  not  singular,  and  the  corresponding 

wave-like  solutions  are  non-singular.  Unless  stated  otherwise,  the  wave-like  solutions  considered 
in  the  following  sections  are  non-singular.  For  these  solutions,  comparing  the  two  results  in  (141) 
generates  the  second  expression  in  (137),  while  the  operations  of  contracting  this  expression  with 
Uj  then  returns  the  first  expression  in  (137),  which  in  view  of  q  =  q(r]i)  can  also  be  expressed  as 


dt 


=  -Hq) 


dq  drji 
dqi  dxe^^ 


(142) 


The  velocity  component  X  =  X{q)  is  determined  by  imposing  the  condition  that  the  non-hnear 
wave-like  expression  (131)  satisfies  the  Euler  equations  (1).  For  q  =  q{r]i),  system  (1)  becomes 


^  dfjdqi  dq 
dt  dq  dxj  drji 


(143) 


Substitution  into  this  system  of  (142)  and  the  second  of  (136),  respectively  for  ||  and  leads 


to 


drji 

dxe 


(144) 


which  yields  the  eigenvalue  problem 


(145) 


B  Wave-Propagation  and  Characteristic  Surfaces 


B.l  Wave-Propagation  Surfaces 

The  wave  propagation  surface  is  the  surface  where  acoustic  and  convection  waves  propagate  and 
therefore,  is  the  surface  that  remains  tangent  to  the  wave-propagation  hyper-planes  expressed  by 
(132)  for  A  corresponding  to  a  wave-propagation  velocity  component.  The  surface 

qe{xjnj  —  Xt)  =  constant  (146) 


for  each  scalar  component  qg  within  q  is  then  itself  tangent  to  a  characteristic  surface  for  an 
appropriate  “constant”,  because 


grad  qi  =  --(ni,n2,n3, -A)-^  parallel  to  (ni,n2,n3, -A) 


(147) 


as  obtained  from  the  second  expression  in  (137). 

Let,  then 

.F(xi,X2,X3,t)  =  constant  (148) 

represent  the  mathematical  function  of  a  wave-propagation  surface.  Prom  vector  analysis  the  vector 
grad  T  remains  perpendicular  to  the  surface  at  each  point.  Since  the  surface  must  be  tangent  to  the 
propagation  plane  (132),  grad  T  must  be  parallel  at  each  point  to  the  unit  vector  (ni,n2,n3,  —A) 
itself  perpendicular  to  plane  (132).  With  rijUj  =  1,  this  condition  yields 


«!  =  T 


M. 

dxi 


jdT  dT  ’  [W 

Y  dxj  dxj  Y 


d:F 

dX2 


dT 

dxj 


=  T 


dr 

dx^ 


A  =  ±- 


dt 


!  dr  dr  ’  j  dr  dr 

Y  dxj  dxj  Y 


(149) 
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Substitution  of  these  expressions  into  the  eigenvalues  (42)  thus  generates  the  separate  propagation- 
surface  differential  equations 

Tt  +  UjJ^x.=0  ,  (Tt  +  UjTx-^  (150) 

where  a  subscript  variable  denotes  differentiation  with  respect  to  that  variable.  The  different  form 
of  these  two  differential  equations  implies  the  existence  of  multiple  propagation  surfaces. 


B.2  Characteristic  Surfaces 


The  surfaces  of  wave  propagation  are  characteristic  surfaces.  By  definition,  a  characteristic  surface 
is  a  surface  across  which  the  solution  q  of  (1)  may  become  discontinuous.  Let  rji  indicate  a 
coordinate  in  the  direction  normal  to  each  facet  of  the  characteristic  surface,  as  shown  in  Figure  15 
In  terms  of  the  partial  derivative  with  respect  to  this  coordinate  direction,  therefore,  system  (1) 
becomes  singular,  or  equivalently,  the  derivative  of  q  in  this  direction  cannot  be  determined  by  the 
system  itself,  where  q  is  not  necessarily  a  wave-like  solution.  By  expressing  the  partial  derivatives 
of  g  in  (1)  in  terms  of  other  coordinates  qj  as 


dq  _  dq  dqi  dq  _  dq  dqi 
dt  dqi  dt  '  dxj  dqt  dxj 

system  (1)  becomes 

(jdqi  dfjjq)  dqi  \  ^  ^  dfjjq)  dqe\  dq 

dt  dq  dxj  j  dqi  ^  y  dt  dq  dxj  j  dqt 


(151) 


(152) 


where  represent  the  cartesian  components  of  a  vector  normal  to  a  facet  of  each 

“.7^  =  constant”  characteristic  surface.  Prom  vector  analysis,  the  components  of  this  normal  vector 
can  then  be  expressed  in  terms  of  the  gradient  of  .7^  as 


dqi  _,dT 
dt  dt 


dqi  _  j^df 
dxj  dxj 


where  A:  is  a  single  scalar  constant.  With  this  specification,  system  (152)  becomes 


k  (i^  +  (i^  +  g/j(g)  dqA 

y  dt  dq  dxj  J  dqi  dt  dq  dxj  J 


dq 

dqe 


(153) 


(154) 


which  becomes  singular,  hence  cannot  determine  when  the  “Ihs”  Jacobian  matrix  is  itself 
singular.  Accordingly,  the  governing  equation  for  the  characteristic  surface  T  is  the  vanishing 
determinant 


dT  ,  ^  dU{q)  dT 


dq  dxj 


=  0 


which  becomes 


Ej)  ^  "I"  UjJ" jj.  ^  C  J" XjJ" I  j  ^  —  0 

This  expression  yields  the  following  partial  differential  equations  for  T 

"t"  ~  0  )  —  C  !F zj^Xj 


(155) 

(156) 

(157) 


identical  to  (150),  which  proves  that  the  smfaces  of  wave  propagation  are  characteristic  surfaces. 
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B.3  Chziracteristic  Cones 

The  principal  propagation  directions  as  well  as  the  supersonic-  and  subsonic-  flow  domain  of  depen¬ 
dence  and  range  of  influence  of  a  flow  field  point  P  in  the  time-space  continuum  can  be  determined 
by  studying  the  variation  of  the  characteristic-surface  shape  versus  the  Mach  number.  To  investi¬ 
gate  these  shape  changes  it  suffices  to  study  the  shape  changes  of  the  corresponding  characteristic 
cones,  which  are  tangent  at  each  flow  field  point  P  to  the  corresponding  characteristic  surfaces, 
but  are  far  easier  to  determine. 

The  characteristic  cones  at  each  flow  field  point  P  are  readily  determined  by  significantly 
simplifying  equations  (157)  through  the  following  Gahlean  space-time  coordinate  transformation 


u^t 

X2  =  X2—  U2t 
Xs  =  Xz  —  Uzt 

T  =  Ct 


(158) 


where  r  now  denotes  a  “space-like”  variable  corresponding  to  the  distance  traveled  by  a  constant- 
speed  aeoustic  wave.  With  this  coordinate  transformation,  the  equation  of  each  eventual  cone 
corresponds  to  the  shape  of  the  cone  as  recorded  by  an  observer  that  moves  along  with  a  fluid 
particle  with  local  convection  velocity  u. 

By  virtue  of  transformation  (158),  the  function  P  is  recast  as 


P(Xi,X2,X3,r)  =P(Xi(xi,t),X2(a;2,t),^3(®3,<),r(t))  (159) 


and  the  partial  derivatives  of  P  with  respect  to  Xi,X2,X3,t  and  xi,X2,  xz,t  are  related  as 


ap 

dP 

BT 

BT 

BP 

BP 

dxi 

X2,X3yt  dXl 

X2,Xz,T  ’ 

Xl,X3,t  QX2 

XuXz,T  ’  dxz 

XlyX2it  dXs 

dt 


Xl,X2,X3 


=  C 


dr 


XuX2,Xs 


dT 

^^dxj 


(160) 


With  these  partial-derivative  relations,  therefore,  the  characteristic-surface  equations  (157)  simplify 
to 


The  equation  of  the  tangent  characteristic  cone  is  then  directly  obtained  from  these  equations 
upon  using  the  correct  multi-dimensional  generalization  of  the  Monge-cone  ordinary  differential 
system,^^"^^  which  yields  the  characteristic-cone  differential  equations 


di^)  d{^)  d{^) 

ag  ,  aG_dz.  ~  M-  j.  dg_^  ~  ,  ag^a^  ~  ag  ,  ag^az  ~ 

axi  ax  axi  ax2  ax  dX2  ax^  ax  axz  ar  ax  at 


dXi  _  dX2  _  dXz  _  dr 

~  __M_  ~  ag  ~  ag  ~  ag 

^ 

The  solution  of  these  ordinary  differential  equations  for  each  equation  in  (161)  is 


(162) 


'  Xi-Xio  =  0 

<  X2-X2o  =  0  ,  iXi-Xio?  +  {X2-X2of  +  {Xz-Xzof  =  {T-Tof  (163) 

:^3-X3o  =  0 
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By  virtue  of  (158),  the  equation  of  the  first  characteristic  cone  with  respect  to  the  fixed  reference 
coordinates  {xi,X2,xz,t)  then  becomes 


xi  -  Xio 

t-to 


=  Ul  , 


X2  -  X2o 
t-to 


=  U2  , 


Xz  —  Xzo 
t-to 


=  Uz 


(164) 


In  the  time-space  continuum  of  a  2-D  flow,  these  equations  correspond  to  the  single  “convection” 
straight  line  S~S'^,  as  displayed  in  Figures  16-17.  As  time  elapses,  therefore,  waves  reach  and  then 
leave  P  along  this  line,  which  projects  onto  the  3-D  streamline 


_  ^2  ^2(?  _  ^3  ^3o 

Ul  U2  Us 

in  the  space  continuum.  A  spatial  discretization  that  aims  to  model  this  propagation  mode, 
therefore,  can  receive  un  upstream  bias  in  the  streamline  direction. 


FigTire  16:  Subsonic  Characteristic  Cone 

In  the  time-space  continuum,  the  equation  of  the  second  characteristic  cone  becomes 
{xi  -  Xio  -Ui{t-  to))^  +  (X2  -  X2o  -  U2{t  -  to))^  +  {X3  -  X30  “  Uz{t  -  <o))^  =  {t  -  tof'  (166) 
hence 

(a;i  -  Xio  -Ui{t-  to)f+{x2  -  X2o  -  U2{t  -  to)f+{x3  -  XSo  -  Uz{t  -  to)f  =  \\uf-j^2  ■ 


(167) 
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This  equation  corresponds  to  a  characteristic  cone  with  vertex  at  P  that  is  tangent  to  the  local 
characteristic  surface  at  P  and  whose  shape  depends  on  M.  As  time  elapses,  therefore,  waves  funnel 
in  towards  P  and  out  and  away  from  P  within  the  characteristic  cone.  In  the  space  continuum, 
expression  (167)  also  represents  the  equation  of  a  family  of  spheres  with  center  coordinates  [xio  + 
Ui{t  —  to)  ,  X2o  +  U2{t  —  to)  ,  xzo  +  uz{t  —  to)]  and  radius  ||u||  {t  —  to)  /M.  For  a  vanishing  c, 
furthermore,  this  cone  collapses  onto  the  convection  line  (164),  which,  however,  is  not  the  axis  of 
the  characteristic  cone.  Figures  16-17  portray  this  characteristic  cone  for  subsonic  and  supersonic 
2-D  flows.  Significantly,  as  M  increases,  the  cone  shears  in  the  streamline  direction,  while  its  cross 
section  remains  an  ellipse,  whose  axes  project  onto  the  streamline  and  crossflow  directions. 

As  a  fundamental  geometric  diflference  between  supersonic  and  subsonic  flows,  a  time  axis 
through  P  is  respectively  outside  and  inside  the  cone.  In  this  manner,  time-independent  rays, 
i.e.  lines  issuing  from  the  time  axis  for  constant  t,  that  are  tangent  to  this  cone  only  exist  when 
the  time  axis  is  not  inside  the  cone,  that  is  for  supersonic  flows.  The  lines  that  issue  from  the 
time  axis  and  axe  tangent  to  spheres  (167)  in  the  {xi,X2,xz)  domain  are  then  found  to  be  Mach 
cones.  This  finding  is  not  coincidental,  because  time-invariant  rays  tangent  to  the  characteristic 
cone  correspond  to  steady-state  ( i.e.  time-  invariant )  characteristic  surfaces,  which  are  the  Mach 
cones  for  the  steady  Euler  equations,  as  remarked  at  the  end  of  Section  4.2. 

Figure  16  presents  the  subsonic  characteristic  cone  for  a  2-D  flow.  The  time  axis  is  inside 
the  cone  and  no  tangent  time-invariant  lines  exist.  Accordingly,  as  time  elapses,  in  this  case,  the 
funneling  of  waves  within  the  cone,  towards  and  then  away  firom  P,  corresponds  on  any  constant-f 
plane  to  closed  curves  encircling  P  that  shrink  towards  P  and  then  expand  away  fi'om  P.  Prom 
this  perspective,  a  region  of  space  centered  about  P  is  simultaneously  the  domain  of  dependence 
and  range  of  influence  for  P.  A  spatial  discretization  that  aims  to  model  this  propagation  mode, 
therefore,  should  receive  un  upstream  bias  along  all  directions  radiating  from  P. 

Figmre  17  presents  the  supersonic  characteristic  cone  for  a  2-D  flow. 
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Figure  17:  Supersonic  Characteristic  Cone 

The  time  axis  is  outside  the  cone  and  the  Mach  lines  are  tangent  to  the  cone  on  each  constant-t 
plane.  As  time  elapses,  in  this  case,  the  tunneling  of  waves  within  the  cone,  towards  and  then 
away  from  P,  corresponds  on  any  constant-t  plane  to  information  reaching  P  within  its  domain  of 
dependence,  between  the  two  Mach  lines  upstream  of  P,  and  leaving  P  within  its  range  of  influence, 
between  the  two  Mach  lines  downstream  of  P.  A  spatial  discretization  that  aims  to  model  this 
propagation  mode,  therefore,  should  receive  an  upstream  bias  along  all  directions  within  the  domain 
of  dependence  of  P,  with  the  streamline  direction  as  the  upstream  bias  principal  direction,  since 
the  streamline  remains  the  axis  of  the  domain  of  dependence  and  region  of  influence. 

The  investigation  of  the  characteristic  cone  shape  change  versus  M,  therefore,  indicates  that 
for  unsteady  and  steady  supersonic  flows,  the  upstream-bias  domain  of  dependence  and  region  of 
influence  for  each  flow  field  point  P  consist  of  the  regions  that  contain  the  streamline  through  P 
and  are  bound  by  the  Mach  lines.  For  unsteady  and  steady  subsonic  flows,  both  the  domain  of 
dependence  and  region  of  influence  encompass  regions  encircling  P. 
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C  Acoustics  Matrices 

Despite  its  zero  eigenvalues,  A^aj  features  a  complete  set  of  eigenvectors  X  and  thus  possesses  the 
similarity  transformation 
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where  the  corresponding  eigenvector  matrix  X  and  its  inverse  X  ^  for  general  equations  of  state 
(7)  axe 
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On  the  basis  of  this  similarity  transformation,  Ajaj  can  then  be  multiply  decomposed  as 

A“  =  XA“+X-^  +  XA“-X-^  ,  A“  =  A“++A“-  (171) 

where  subscript  n  here  denotes  the  number  of  non-zero  eigenvalues  within  An"*”  and  An” .  At  least 
two  specific  decompositions  of  AjOj  can  be  considered,  which  respectively  correspond  to  n  =  1 
and  n  =  3  and  both  lead  to  a  conveniently  simple  bi-modal  upstream  approximation  of  AjOj.  For 
n  =  1,  the  expressions  for  An'*’  and  An“  are: 
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the  corresponding  expressions  for  n  =  3  are 
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For  n  =  3,  the  matrix  product  of  lA^aj  with  the  directional  derivative  Okdq/dxk  of  g  directly 
leads  to  the  beautifully  simple  acoustic-field  result 


U?o,-  =  XcIX-\k-^  =  car 
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The  similar  results  for  (af ^ , 03 ^ )  and  02^,03^)  replacing  (01,02,03)  are 


^2  N2  N2. 


3  *  ^  dXk~  ^  dXk 


yl?of2  O^- 

3  “fc  ^ 


'-1„JV2_^  _  „JV2 


X^cIX-^aZ^i^  =  caZ-^i^ 

^2  axfc  *=  arcfc 


where  /  denotes  the  5x5  identity  matrix.  For  a  vanishing  M,  therefore,  all  the  eigenvalues  of  the 
streamline  and  crossflow  absolute  acoustic  matrices  approach  +c. 


D  Upstream-Bias  Stability  Condition 

In  Jacobian  form,  expression  (88)  becomes 

S = ^  ^  ^  ^ "  “'S "  S 

By  virtue  of  a  non-singular  (a:i,X2,a:3)-space  coordinate  transformation 

rii  =  XjUj  -  X{q)t  ,  m  =  V2ixii^2,X3,t)  ,  m  =  T)3{xi,X2,X3,t)  (178) 

analogous  to  (133),  and  the  inverse  {J‘  =  {^jk}  /  det  J  of  the  associated  coordinate  -  trans¬ 

formation  Jacobian  J  =  dx/dr],  expression  (177)  for  each  t  becomes^^ 

9fj  1  ^  [  /  f  f  ,  jv  z'  jv,  jv,  ,  w,  nA\  t  ,  ^jk  dq 

“  ^‘diTjaS  r‘“  (' “  (“<  +  “i  “»  ))  ^ + “<  a^-  +  j 

(179) ' 

where  the  coordinate-transformation  metric  data  Cjjt?  within  J,  can  be  brought  inside  the  partial 
differential  operation  because  of  the  fundamental  equality^^  dea/drn  =  0,  for  any  i  and  with 
implied  summation  on  t  For  each  moreover,  the  metric  data  ea  denote  the  components  of  a 
vector  that  is  orthogonal  to  an  %  =  constant  surface. 

As  delineated  in  Appendix  A,  wave-like  solutions  q  depend  upon  7]i  only.  For  these  solutions, 
therefore,  q  =  q{r)i)  and  (179)  exhibits  the  upstream-bias  expression 

("  ^ dfj^ 

(180) 

for  which  en / ^/ {eneei)  denotes  the  component  rii  of  the  wave-propagation  unit  vector  n.  This 
expression  now  corresponds  to  a  one-dimensional  upstream-bias  expression.  This  expression  will, 
therefore,  be  stable  when  all  the  eigenvalues  of  the  upstream-bias  matrix 

A^ni(c  (aaittj  +  I  -t-  o-i-^  +  (1®^) 


remain  positive  for  all  M  and  propagation  directions  n. 
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E  Crossflow  Upstream  Function 


For  the  determination  of  =  a^{M)  from  A3^4  (M,n)  >  0  in  (97),  heed  that  a  =  v  remains 
perpendicular  to  both  and  a^2.  Moreover,  the  expressions  njOi,  njofS  and  respectively 
denote  the  vector  “dot”  products  between  the  unit  vector  n  and  the  unit  vectors  a,  ,  and  ^ 
hence 

n  =  njOja  +  (182) 

Accordingly, 

riitti  =  cos  0 

(uiaf^aj^rij  +  nia^^aj^nj^  =  |^n  •  = 

=  |P  cos^(90  —  6)  nia^^aj^rij  +  Uiaf'^a^^nj  =  sin^  6 

where  6  denotes  the  angle  between  n  and  a. 

For  eigenvalue  A4  =  A4  {M,n)  in  (94),  therefore,  the  condition  A4  (M,n)  >  0  yields 


a^>9{0,M)  = 


cos  0y^^-^-^PgCos0M^  4^  —  cos^  0  +  ^1  +  ~^Pe^ 


1  —  cos2  6 

For  supersonic  flows  with  M  >  1  +  a  =  0  and  <S  =  1,  hence  (184)  becomes 

cos9  —  M  cos^O 


a">9(e.M)  = 


(183) 


1  —  cos^  9 


(184) 


(185) 


and  in  particular  >  5max(Af),  where  ^tnax(Af)  denotes  the  maximum  oi  g  =  g{6,M)  with 
respect  to  6,  for  each  M.  From  (185),  the  determination  of  gmax{M)  yields 

^  =  0  =»  cos^  6  —  2M  cos  0  +  1  =  0  (186) 

which  leads  to 

cos^l  =  M  -  VM2  -  1  ,  5max(M)  =  I  (m  -  \/M2  -  l)  (187) 

•9 — 5'inax  ii  '  ' 

Significantly,  the  same  solution  for  gmax{M)  results  from  the  condition  A3  (M,n)  >  0.  Conse¬ 
quently, 

a’^  >^{M-  \/M2  -  1)  ,  M>Mm  =  H-£m  (188) 

and  considering  that  A4  (1,  w)  =  an  analogous  equality  is  adopted  for  leading  to 

=  Qmaxi^M)  H"  ^m/2-  ^ 

For  subsonic  flows,  a  numerical  analysis  of  3  =  g{6,  M)  from  (185)  reveals  that  g{0^  M)  <  0.3  for 
all  0  and  M  <  Mm-  Additionally,  for  an  isotropic  acoustic  upstream  for  vanishing  M,  (0)  =  1 
and  da^ /dM\M=o  =  0,  whereas  for  M  >  Mm  both  and  its  derivative  with  respect  to  M 
follow  from  (188).  A  smooth  variation  for  =  a^{M)  that  satisfles  all  of  these  constraints  is  the 
composite  spline 
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where  superscript  prime  “ '  ”  denotes  differentiation  with  respect  to  M. 
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(189) 
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F  Discrete  Conservation 


For  a  clear  comparison  between  finite  difference/ volume  schemes  and  the  characteristics-bias  finite 
element  algorithm  (109),  the  specific  expressions  for  the  several  averages  in  (115)-(116)  are 
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for  the  flux  averages, 
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for  the  average  streamwise  acoustics  contributions  to  Fi  .  and  F2  , .  The  corresponding  cross- 

»+5J  ‘’^  +  f 

flow  acoustics  contributions  directly  follow  from  these  expressions  by  replacing  (a,ai,a2)  with 
,02  ).  The  flux  divergence  averages  are  then  expressed  as 
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and 
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which  provide  the  contributions  to  Fi.  ,  and  F^. .  ,  from  the  average  divergence  of  the  flux  /?. 

The  corresponding  contributions  from  the  divergence  of  the  flux  are  then  directly  obtained  from 
these  expressions  by  replacing  with  ff  and  (01,02)  with  (oi(5, 02^).  In  all  these  averages,  the 
ratios  (£)e/Axi  and  (£)e/Ax2  are  positive  numbers  of  order  1,  for  (^)e  remains  commensurate  with 
\/(Axi)2  +  (Ax2)^. 

Equation  (113)  along  with  expressions  (115)-  (116),  (190)-(195)  can  also  be  rearranged  as  a 
linear  combination  of  two-point  upstream  and  downstream  flux  differences.  For  instance,  the  linear 
combination  L^-j  of  the  differences  of  along  the  segment  i  —  l,i,i  -\-  i  on  level  j  of  elements 
i  -  ^,j  -  5  and  i  -t-  5,  j  -  5  in  Figure  11,  becomes 
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In  this  and  the  analogous  linear  combinations  of  flux  differences,  the  numerical  values  of  the  con¬ 
troller  determine  the  combination  weights  of  the  downstream  and  upstream  expressions.  Since 
remains  non-negative,  these  equations  induce  an  upstream  bias  on  the  appropriate  flux  differ¬ 
ences,  for  the  weight  on  the  upstream  flux  difference  is  greater  that  the  weight  on  the  downstream 
difference.  For  example,  for  02  =  0  and  oi  =  1  then  {£)e  =  Axi  and  the  weight  on  the  downstream 
flux  difference  vanishes  altogether  for  =  !•  As  a  result,  the  finite  element  weak  statement  (112) 
generates  consistent  variable-upstream-bias  discrete  equations  that  correspond  to  an  upstream-bias 
discretization  of  the  original  Euler  system  (1). 
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SUMMARY 

This  paper  introduces  a  continuum,  i.e.  non-discrete,  upstream-bias  formulation  that  rests  on  the  physics 
and  mathematics  of  acoustics  and  convection.  The  formulation  induces  the  upstream  bias  at  the  differential- 
equation  level,  within  a  characteristics-bias  system  associated  with  the  Euler  equations  with  general  equi¬ 
librium  equations  of  state.  For  low  subsonic  Mach  numbers,  this  formulation  returns  a  consistent  upstream- 
bias  approximation  for  the  non-linear  acoustics  equations.  For  supersonic  Mach  numbers,  the  formulation 
smoothly  becomes  an  upstream-bias  approximation  of  the  entire  Euler  flux.  With  the  objective  of  mini¬ 
mizing  induced  artificial  diffusion,  the  formulation  non-linearly  induces  upstream-bias  essentially  locally  in 
regions  of  solution  discontinuities,  whereas  it  decreases  the  upstream-bias  in  regions  of  solution  smoothness. 
The  discrete  equations  originate  from  a  finite  element  discretization  of  the  characteristic-bias  system  and  are 
integrated  in  time  within  a  compact  block  tri-diagonal  matrix  statement  by  way  of  an  implicit  non-linearly 
stable  Runge-Kutta  algorithm  for  stiff  systems.  As  documented  by  several  computational  results  that  re¬ 
flect  available  exact  solutions,  the  acoustics-convection  solver  induces  low  artificial  diffusion  and  generates 
essentially  non-oscillatory  solutions  that  automatically  preserve  a  constant  enthalpy  as  well  as  smoothness 
of  both  enthalpy  and  mass  flux  across  normal  shocks. 

No,  of  Figures:  24.  No.  of  Tables:  1.  No.  of  References:  13. 

KEY  WORDS:  CFD;  upwind;  artificial  diffusion;  acoustics;  convection;  finite  elements;  implicit  integration 

1  Introduction 

This  paper  introduces  for  the  Euler  equations,  with  general  equilibrium  equations  of  state,  a  stable 
and  consistent  upstream-bias  algorithm  that  rests  on  a  new  Flux  Jacobian  Decomposition  (  FJD  ) 
formulation,  features  the  simplicity  of  a  flux  vector  splitting  formulation  and  accommodates  an 
implicit  solver.  The  algorithm  induces  minimal  diffusion,  naturally  incorporates  a  finite  element 
discretization  and  uniquely  generates  the  upstream  bias  directly  at  the  difierential  equation  level 
before  and  independently  of  any  discrete  approximation  on  specified  grids.  As  one  of  its  impor¬ 
tant  features,  the  algorithm  combines  the  mathematics  of  upwind  algorithms  with  the  physics  of 
acoustics  and  convection,  the  wave  propagation  mechanisms  within  gas  dynamic  flows. 

Most  finite  element,  difference,  and  volume  algorithms  have  remained  largely  independent  fi:om 
the  physics  of  acoustics  and  convection.  The  dissipation  mechanisms  within  these  algorithms, 
furthermore,  have  been  developed  at  the  discrete  level  in  connection  with  a  specific  grid. 

Several  finite  element  Euler  solvers  have  either  utilized  modifications  of  the  test  space,  e.g. 
SUPG^,  or  introduced  Taylor’s  series  based  dissipation  terms,  e.g.  TWS^,  to  generate  stable 
algorithms.  The  mathematical  developments  in  these  fundamental  contributions  have  remained 
independent  from  characteristics  theory. 

Most  characteristics  upwind  schemes  have  been  implemented  through  finite  difference  and  finite 
volume  discrete  approximations^.  Roe’s  approximate  Riemann  solver"^  remains  an  upwind  proce¬ 
dure  with  a  significant  physics  content.  The  associated  operation  count,  on  the  other  hand,  exceeds 
that  of  a  flux  vector  splitting  formulation. 

*  Associate  Professor,  jicinnell@utk.edu.  Joe  lannelli,  Dept,  of  Mechanical  &  Aerospace  Engineering  and  Engi¬ 
neering  Science,  The  University  of  Tennessee,  315  Perkins  Hall,  Knoxville,  TN  37996-2030,  U.S.A. 


2 


Joe  lannelli 


VanLeer’s  original  simple  flux  vector  splitting  and  the  many  variants  developed  thereafter® 
essentially  rely  on  Mach-number  dependent  polynomials  to  generate  flux  components,  each  featur¬ 
ing  Jacobian  eigenvalues  with  uniform  algebraic  signs.  As  such,  these  fundamental  mathematical 
developments  remain  independent  from  the  physics  of  acoustics  and  convection. 

As  a  spin-off  from  these  studies  and  presented  as  a  new  perfect-gas  flux  vector  splitting,  Liou’s 
and  Steffen’s  procedure®  employs  an  ad-hoc  advection  velocity  along  with  flux  components  with 
convection  and  pressure  physical  meanings.  An  analogous  method  was  introduced  earlier  by  this 
author^  and  the  essentially  non-oscillatory  results  from  both  methods  bear  out  the  advantages  of 
employing  flux  components  with  clear  physical  significance. 

The  eigenvalues  of  the  various  convection  and  pressure  fluxes  in  these  methods,  however,  either 
approach  zero  or  remain  substantially  less  than  the  speed  of  sound  for  decreasing  Mach  number. 
Neither  method  has  therefore  generated  a  physically  consistent  upstream-bias  approximation  of  the 
acoustics  limit  of  the  Euler  equations  in  the  low-Mach  number  regime.  No  single  decomposition  of 
the  Euler  flux,  in  fact,  contains  separate  components  that  respectively  correspond  to  the  physics 
of  acoustics  and  convection. 

In  this  paper,  a  new  upstream-bias  formulation  is  developed  that  is  based  on  a  decomposition  of 
the  Euler  flux  vector  Jacobian  into  matrix  components.  This  formulation  encompasses,  unifies  and 
generalizes  upwind  algorithms,  including  Flux  Vector  Splitting  and  Flux  Difference  Splitting  devel¬ 
opments.  This  formulation  develops  the  upstream-bias  approximation  directly  at  the  differential 
equation  level,  before  any  discretization.  The  method  results  in  a  “companion”  characteristics- 
bias  system  that  is  associated  with  the  Euler  equations  and  contains  an  upstream-bias  differential 
expression.  The  characteristics-bias  system  features  a  characteristics  flux  that  generalizes  in  the 
continuum  the  traditional  numerical  fluxes  of  upwind  schemes.  The  acoustic-convection  upstream 
algorithm  then  results  from  a  specific  decomposition  of  the  flux  vector  Jacobian  into  genuine  acous¬ 
tics  and  convection  components,  for  a  physically  consistent  upstream  approximation  of  coupled 
acoustic  and  convection  wave  propagation. 

A  traditional  centered  discretization  of  the  acoustics-convection  characteristics-bias  system  then 
automatically  generates  a  coherent  upstream  discrete  approximation  of  the  governing  Euler  equa¬ 
tions.  This  approximation,  moreover,  reduces  to  a  consistent  upstream  approximation  of  the  acous¬ 
tics  equations,  for  vanishing  Mach  number,  which  addresses  the  challenging  problem  of  calculating 
low-Mach-number  flows.  Finite  difference,  volume,  or  element  procedures  can  be  used  to  discretize 
the  characteristics-bias  system.  The  algorithm  in  this  paper  has  used  a  finite  element  discretiza¬ 
tion,  which  also  leads  naturally  and  automatically  to  consistent  boundary  differential  equations 
and  a  new  outlet  pressme  boundary  condition  that  does  not  require  any  algebraic  extrapolation  of 
variables.  The  resulting  discrete  equations  correspond  to  an  essentially  centered  discretization  in 
the  form  of  a  non-linear  combination  of  upstream  diffusive  and  downstream  anti-diffusive  flux  dif¬ 
ferences,  with  greater  bias  on  the  upstream  diffusive  flux  difference.  This  formulation,  furthermore, 
directly  accommodates  an  implicit  solver,  for  it  is  very  easy  to  determine  the  required  Jacobian 
matrices. 

With  the  objective  of  minimizing  induced  artificial  diffusion,  the  characteristics-bias  flux  non- 
linearly  induces  upstream-bias  essentially  locally  in  regions  of  solution  discontinuities,  whereas  it 
decreases  the  upstream-bias  in  regions  of  solution  smoothness.  This  variable  diffusion  is  auto¬ 
matically  adjusted  at  each  discretization  node  by  a  new  controller,  also  introduced  in  this  paper. 
This  controller  depends  on  local  solution-slope  Jumps  and  varies  the  combination  weights  on  the 
upstream  and  downstream  fluxes,  within  the  discrete  equations. 

The  operation  count  for  this  algorithm  is  comparable  to  that  of  a  simple  flux  vector  splitting 
algorithm.  The  developments  in  this  study  have  employed  basic  two-noded  cells,  which  has  thus 
led  to  a  block  tridiagonal  matrix  system,  for  the  implicit  formulation.  To  determine  the  ultimate 
accuracy  of  linear  approximations  of  fluxes  within  two-noded  cells,  for  a  computationally  efficient 
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implementation,  this  study  employs  no  MUSCL-type  local  extrapolation  of  dependent  variables. 

This  paper  is  organized  in  9  sections.  After  the  introductory  remarks  in  Section  1,  Section  2 
presents  the  governing  equations  and  Sections  3  develops  the  Flux  Jacobian  Decomposition  for¬ 
mulation.  Section  4  presents  the  decomposition  of  the  Euler  flux  Jacobian  into  genuinely  physical 
acoustic  and  convection  components,  followed  in  Section  5  by  the  determination  of  the  correspond¬ 
ing  upstream-bias  stability  eigenvalues.  Section  6  details  the  finite  element  spatial  discretization, 
along  with  the  new  pressure  boundary  condition  and  solution  dependent  control  of  upstream  diffu¬ 
sion,  and  Section  7  delineates  the  non-linearly  stable  implicit  Runge-Kutta  time  integration.  The 
computational  results  are  discussed  in  Section  8,  with  concluding  remarks  presented  in  Section  9. 

2  Governing  Equations 

2.1  Euler  System 

With  respect  to  an  inertial  reference  frame,  the  quasi-lD  Euler  conservation  law  system^  is: 


9q  df{q) 


(1) 


where  the  independent  variable  {x,t)  varies  in  the  domain  D  =  Q,x  [to)T],  =  [a,  6].  This  system 

consists  of  the  continuity,  momentum  and  total-energy  equations,  and  the  arrays  q  =  q{x,t), 
f  =  f{q)  and  ^  =  ^{x,q)  are  defined  as 
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where  p,  m,  E  respectively  denote  static  density  and  volume-specific  linear  momentum  and  total 
energy;  the  eulerian  flow  velocity  u  is  then  defined  as  u  =  m/p.  The  terms  p,  A  =  A{x)  and  A* 
respectively  indicate  static  pressure,  area  of  the  flow  duct  cross-section,  and  constant  upstream-flow 
reference  throat  area. 


2.2  Equilibrium  Equation  of  State,  Pressure  Derivatives,  Speed  of  Sound 

For  any  homogeneous  equilibrium  gas,  pressure  depends  upon  two  other  thermodynamic  variables.® 
They  are  density  p  and  mass-specific  internal  energy  e,  in  this  case,  since  they  are  readily  available 
from  the  Euler  system  (1):  p  directly  from  the  continuity  equation  in  the  system,  and  e  from  q  as 


=  E 

~  P 


2p‘ 


rm 


(3) 


(4) 


The  pressure  equation  of  state  thus  becomes 

The  Jacobian  derivatives  of  p  with  respect  to  g,  for  the  Jacobian  df/dq  of  f{q),  are  not  all  in¬ 
dependent  of  one  another.  The  derivatives  of  (4)  with  respect  to  m  and  E  in  fact  satisfy  the 
constraint 


dp 
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as  obtained  by  expressing  the  derivatives  of  p  with  respect  to  m  and  E  in  terms  of  the  thermody¬ 
namic  derivative  of  p  with  respect  to  e,  from  the  first  expression  in  (4).  In  the  following  sections, 
for  simplicity,  the  abridged  notation 


dp 


Pm  = 


dp 

dm 


p,E 


Pe 


dp 


P,TTl 


(6) 


will  denote  the  Jacobian  derivatives  of  pressure.  The  particular  perfect-gas  expressions  for  (4) 
follow  from  the  internal  energy  and  equation  of  state 


e  =  cr=-^T  ,  p  =  pRT  (7) 

7-1 

for  this  type  of  gas,  where  c^,,  T,  and  7  respectively  denote  the  constant-volume  specific  heat, 
static  temperature,  gas  constant  and  specific-heat  ratio.  The  elimination  of  T  from  these  two 
expressions  and  use  of  (3)  leads  to  the  following  familiar  expressions  for  the  equation  of  state  for 

p  in  terms  of  q  /  1  \ 

p  =  {'r-l)pe  =  {7-l)i^E-—m^j  (8) 

The  square  of  the  speed  of  sound  c  for  general  equilibrium  equations  of  state  can  be  expressed 


as 


=  Pp 

s 


(9) 


in  terms  of  the  Jacobian  partial  derivatives  of  p.  With  this  result,  the  mass-specific  total  enthalpy 
H  depends  on  g  as 


(10) 


where  M  =  ||m||/c  denotes  the  Mach  number. 


2.3  Characteristics  Analysis 

For  general  equilibrium  pressure  equations  of  state  (4),  the  characteristic  speeds  associated  with 
the  Euler  equations,  that  is  the  eigenvalues  of  the  flux  vector  Jacobian 
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have  been  exactly  determined  in  closed  form  as 


Ai  =  u 


A2,3  =  u  ± 


1/2 


(11) 


(12) 


These  eigenvalues  correspond  to  the  slopes  of  the  characteristics,  as  portrayed  in  Figure  1  for 
representative  hypersonic,  supersonic,  sonic,  and  subsonic 

Of  interest,  eigenvalues  A2,3  directly  incorporate  a  sound  speed  expression  that  coincides  with 
the  isentropic  partial  derivative  of  pressure  (9).  Through  (9),  therefore,  these  equilibrium-gas 
eigenvalues  become 

-  u  ,  A?,  =  u±c  (13) 
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Figure  1:  Characteristics:  a)  Hypersonic,  b)  Supersonic,  c)  Sonic,  d)  Subsonic 


which  have  the  same  familiar  form  as  the  perfect-gas  eigenvalues. 

Figure  1  shows  the  characteristics  in  a  suitable  neighborhood  of  a  flow  field  point  P  in  a  (t,  x) 
plane.  An  interesting  geometric  diff“erence  among  supersonic,  sonic,  and  subsonic  flows  is  that  a 
time  axis  through  P  is  respectively  outside,  on  the  boundary,  and  inside  the  domain  of  dependence 
and  range  of  influence  of  point  P.  Wave  propagation  for  supersonic  flows  essentially  occurs  by 
convection,  mono-axially  from  upstream  to  downstream  of  P;  the  sonic  case  becomes  a  limiting  case; 
for  subsonic  flows,  instead,  wave  propagation  occurs  by  both  convection  and  acoustics,  bi-modally 
from  both  upstream  and  downstream  toward  P;  for  vanishing  Mach  number,  wave  propagation  is 
essentially  acoustic. 

Since  gas  dynamic  wave  propagation  physically  occurs  by  acoustics  and  convection,  the  up¬ 
stream  CFD  algorithm  in  this  paper  mathematically  models  this  coupled  acoustic- convection  wave 
propagation.  The  algorithm  identifies  the  genuine  convection  and  acoustics  components  within  the 
flux  Jacobian  and  then  establishes  a  physically  consistent  upstream  approximation  for  each  of  these 
components. 

2.4  Non-Linear  Acoustics  Equations 

The  Euler  equations  contain  the  acoustics  equations  for  vanishing  Mach  numbers.  Identification 
of  these  equations  yields  the  acoustics  component  of  the  Euler  flux  Jacobian  for  any  Mach  number. 
Upon  writing  momentum  m  in  terms  of  the  Mach  number  M  as  m  =  pcMu/\u\  and  using  the 


energy  pressure  derivative  identity  (5),  the  Euler  system  (1)  becomes 


1 


dt  dx  u|  A  dx 


,  do  dE  u  \,^  .dm  mdp  mdA'\  _ 
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dE  E  +  pdm  ^  \  d  (EAp^  EAy^ 


A  dx 


and  for  a  vanishing  Mach  number,  these  equations  reduce  to  the  non-linear  acoustics  system 


P  \ 

m  + 
E  } 


^  ^  ^  S  ^ 

0  )  Pe  -  771  I  =  0 

E  +  P  .  5a;  p 


where  A"’  denotes  the  rero-Mach-number  aeoustics  matrix.  The  dependent  variables  in  th.^  equa, 
tacmrespL  to  those  in  a  flow  field  that  originates  from  slight  perturbafons  to  an  otherw.se 

‘*“S‘ttuhe  energy  equation  toward  steady  state,  in  this  case,  is  no  longer  linearly  independent 
fromrhe  “inuity  equation.  This  phenomenon  directly  explains  the  widely  reported  convergent 
drculfa  experieLd  in  the  CFD  simulation  of  incompressible,  i.e.  low-Mach-number,  flows  w,th 

a  compressible  flow  formulation.  .an 

By  virtue  of  the  total  enthalpy  expression  (10),  the  matrix  A  becomes 


A“<>  = 


(T  -P, 


with  eigenvalues  ,  n?! 

A“'>  =  0  ,  A“°3  =  ±c 

mere  c  corresponds  to  the  aero-Mach-number  isentropic  speed  of  sound  (9).  With  XT  =  0,  the 
propagation  of  (acoustic)  waves  governed  by  this  system,  therefore,  corresponds  to  an  isentrop 

process  with  negligible  flow  kinetic  energy. 

Equations  (15)  contain  the  important  expressions 

dE  _  -  Pf,  dp  dp  Pe  dE 

~  Pg  dx  '  dxc^  -  Pp  dx 

which  result  from  the  momentum  equation  by  expressing  p,  therein  using  (9)  as 


dp  dp 

~dx~  ^^dx 


fdp  E  +  p\dp  f— -  (19) 

~\fps~^^  P  j  5a:  dx  dp  s  dx  ^^\dx  p  dx) 

For  an  isentropic  flow  p  =  p(p),  hence  the  first  rhs  term  in  (19)  equals  |f .  The  second  rhs  term  must 
consequently  vanish,  which  returns  results  (18),  after  using  (10).  These  results  will  convenien  y 
simplify  the  acoustics-convection  upstream  formulation  in  Sections  4.2-  .  . 
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3  Non-Discrete  Upstream-Bias  Approximation 


The  non-discrete,  i.e.  continuum  or  before  discretization,  upstream-bias  approximation  of  the  Euler 
equations  derives  from  a  characteristics-bias  integral  statement  associated  with  (1).  The  prototype 
integral  statement  is 

r  /  fi/y 

(20) 


.dt 

which  is  equivalent  to  the  governing  system  (1)  for  arbitrary  subdomains  C  fi  and  arbitrary  test 
functions  w  with  compact  support  in  The  characteristic-bias  integral  is  then  defined  as 


L 


^  ( dq  df^' 


da  =  o 


(21) 


where  corresponds  to  a  characteristics-bias  flux  that  automatically  induces  within  (21)  an 
upstream-bias  approximation  for  the  Euler  flux  divergence 

3.1  Flux  Jacobian  Decomposition  and  Upstream-Bias  Integral  Average 
To  develop  the  flux  /^,  consider  first  the  flux  Jacobian  decomposition  (  FJD  )  into  L  contributions 


^  A 


e=i 


di 

dx 


=  ^  at  At 


l=\ 


dx 


(22) 


where  At  corresponds  to  a  flux-jacobian  matrix  component  with  uniform-sign  eigenvalues  and  at 
denotes  a  linear-combination  function,  possibly  depending  upon  q. 

An  integral  average  of  the  Euler  flux  divergence  ^  as  expressed  through  decomposition  (22) 
becomes 

(23) 


Lw—d9,=  ('^watAt^d'^ 

Jn  ox  Jn  ^  ox 

The  flux  is  therefore  defined  by  way  of  an  upstream-bias  integral  average  as 


(24) 


where  the  rhs  integral  provides  an  upstream  bias  for  each  matrix  component  within  the  FJD  in 

(22). 

The  positive  'ip  in  (24),  0  <  ■0  <  1,  stands  for  a  new  “upstream-bias”  controller,  which  au¬ 
tomatically  adjusts  the  amount  of  induced  upstream-bias  diffusion,  depending  on  local  solution 
non-smoothness,  as  introduced  and  detailed  in  Section  6.3.  The  variation  5tw  induces  the  appro¬ 
priate  upstream-bias  for  the  test  function  w  for  each  component  within  (24).  Depending  on 
the  physical  significance,  magnitude  and  algebraic  sign  of  the  eigenvalues  of  At,  the  variation  6tw 
can  vanish  or  become  algebraically  positive  or  negative,  which  corresponds  to  an  upstream  bias 
respectively  in  the  negative  or  positive  sense  of  the  x-axis. 


3.2  Characteristics-Bias  Flux 
The  variation  5tw  in  (24)  becomes 

.  ^  dw .  dw 

Otw  =  -^otx  =  -^ate  ,  ate  =  Otx 
dx  dx 


(25) 
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where  e  denotes  a  local  positive  length  scale,  while  the  direction  cosine  ai  can  equal  0  or  +1,  —  1, 
possibly  also  depending  upon  q. 

With  these  specifications,  the  upstream-bias  integral  average  (24)  becomes 


/- 


dx 


(26) 


Considering  that  w  has  compact  support  in  fi,  it  vanishes  on  the  boundary  dQ  of  fl.  As  a  result, 
integrating  (26)  by  parts  generates 


dx  dx 


\  fci 


=  0 


(27) 


which  contains  no  boundary  integrals.  Since  this  integral  must  vanish  for  arbitrary  test  functions  w 
and  domains  fi,  its  integrand  must  identically  equal  zero,  which  generates  the  following  expression 
for  the  divergence  of  the  chaxacteristics-bias  fiux 


df  ^  df 

dx  dx 


dx 


I  e'lp'^aeatAi 
\  i=i 


(28) 


This  expression  exhibits  an  upstream-bias  artificial  difiPusion,  in  the  form  of  a  second-order  differ¬ 
ential  expression  with  matrix 

L 

^  ^  agaeAe  (29) 

^=1 

For  physical  consistency  of  the  upstream  bias  in  (24)- (28)  and  associated  mathematical  stability  of 
the  corresponding  second-order  differential  expression,  all  the  eigenvalues  of  this  upstream  matrix 
must  be  positive.  This  requirement  becomes  a  fundamental  upstream-bias  stability  condition. 

For  2-D  and  3-D  flows,  the  version  of  the  characteristics-bias  flux  follows  as  a  multi-dimensional 
generalization  of  (22)- (29)  as 


and 


dfj  d  (  ,4^  ,  dq\ 


(30) 


(31) 


with  1  <  i,j  <  2,  for  2-D,  and  1  <  i,j  <  S,  for  3-D.  In  these  expressions,  an  now  denotes  the 
t***  direction  cosine  of  a  unit  vector  along  the  principal  wave-propagation  direction  of  wave  ’ , 
and  LUi  indicates  the  direction  cosine  of  a  unit  vector  u>  along  an  arbitrary  wave-propagation 
direction.  The  multi-dimensional  acoustics-convection  form  of  (31)  is  being  completed,  but  will 
be  detailed  soon  in  a  companion  paper,  to  keep  the  length  of  each  of  these  two  papers  within  an 
appropriate  size. 

The  continuum  expression  (28),  or  (31),  for  the  divergence  of  the  characteristics-bias  flux  consti¬ 
tutes  a  non-discrete  generalization  of  the  various  numerical-flux  formulae  employed  in  several  CFD 
upwind  schemes.  It  encompasses,  generalizes,  and  unifies  flux-vector  and  flux-difference  schemes 
as  shown  by  the  following  representative  examples. 


Acoustics-Convection  Euler  Solver 


9 


3.2.1  van  Leer’s  Formulation  and  Flux  Vector  Splitting 

Consider  the  van  Leer’s  formulation  as  a  representative  Flux  Vector  Splitting  (  FVS  ).  In  this 
formulation,  the  inviscid  flux  /  is  “split”  as 

/  =  (32) 

where  the  Jacobian  matrices  of  and  respectively  possess  non-negative  and  non-positive 
eigenvalues. 

The  FJD  expression  (22)  encompasses  (32)  with  L  =  2  as 


L  f)fVL+  FifVL- 

y)  oieAt  =  — ^ - 1 - ^ —  ,  Q!i  =  1  ,  02  =  1  (33) 

The  corresponding  characteristics-bias  flux  divergence  for  van  Leer’s  FVS  accrues  from  (28) 
with  i/}  =  I,  ai  =  1,  02  =  —1  as 

f  !?)  =  £/  9  ,34, 

ox  ox  ox  \  y  oq  oq  I  ox  J  ox  ox  y  y  ox  ox  j  J 


which  generalizes  in  the  continuum  the  traditional  numerical  flux  formulae  for  FVS  constructions. 
The  associated  upstream  matrix  A  is 


gyVL+  QfVL- 

dq  dq 


(35) 


The  upstream-bias  stability  condition,  however,  is  not  automatically  satisfied,  even  though  each  of 
the  two  matrices  ^  and  ^  has  positive  eigenvalues.  This  stability  condition  is  not 

unconditionally  satisfied  because  the  sum  of  two  positive-eigenvalue  matrices  does  not  necessarily 
yield  a  matrix  with  positive  eigenvalues.  As  an  example  consider  the  following  matrix  sum  of  two 
positive-eigenvalue  matrices 


where  a  is  a  real  number.  One  of  the  eigenvalues  of  this  matrix  sum  is  negative  for  cr  >  4.  For 
instance,  for  a  =  7  the  eigenvalues  are  4-9  and  —1. 

Most  likely,  however,  (35)  satisfies  the  upstream-bias  stability  condition  for  most  of  the  flow 
conditions  considered  in  the  technical  literature,  in  view  of  the  stable  results  reported.  For  sub¬ 
sonic  flows,  each  of  the  two  flux  vector  components  in  (32)  remains  unrelated  to  the  physics  of 
acoustics  or  convection.  On  the  other  hand,  (32)  is  computationally  advantageous,  for  it  calls  for 
the  discretization  of  simple  flux-vector  components. 


3.2.2  Roe’s  Formulation  and  Flux  Difference  Splitting 

Consider  next  Roe’s  formulation  as  a  representative  Flux  Difference  Splitting  (  FDS  )  development. 
In  this  formulation,  the  inviscid  flux  Jacobian  of  /  is  “split”  as 


^  =XA+X-^  -[-XA-X-^ 

dq 


(37) 
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where  X  and  A  =  A+d-A"  denote  the  right  eigenvector  matrix  and  eigenvalue  diagonal  matrix  of  the 
iacobian,  all  evaluated  at  special  average  values  of  g,  with  A+  and  A  respectively  containing  non¬ 
negative  and  non-positive  eigenvalues.  The  matrices  at  the  rhs  of  (37),  therefore,  will  respectively 
possess  non-negative  and  non-positive  eigenvalues. 

The  FJD  expression  (22)  encompasses  (37)  with  L  =  2  as 


=  XA+X"^  +  XA  X  ^  ,  ai  =  1  ,  0:2  =  1  (38) 

The  corresponding  characteristics-bias  divergence  for  Roe’s  formulation  accrues  from  (28)  with 

^  =  1,  ai  =  1,  02  =  “"1  ^ 

^  =  I  -  S  ('  s )  =  s  -  i 

which  generalizes  in  the  continuum  the  traditional  numerical  flux  formulae  for  FDS  constructions. 
The  associated  uspstream  matrix  A  is 

^  =  X(A+-A-)X-'  (40) 

which  has  non-negative  eigenvalues  and  therefore  automatically  satisfies  the  upstream-bi^  stability 
condition  for  any  flow  state  for  which  no  eigenvalue  vanishes.  The  discretization  of  (39)  calls  tor 
more  computational  operations  than  (34),  while  each  of  the  two  rhs  components  in  (37)  lumps  into 
one  matrix  the  matrices  representative  of  the  distinct  acoustics  and  convection  wave  propagation 
mechanisms.  On  the  other  hand,  numerous  numerical  results  bear  out  the  accuracy  of  an  h  Ub 

formulation. 


4  Acoustics-Convection  Flux  Jacobian  Decomposition 

The  acoustics-convection  flux  Jacobian  decomposition  consists  of  components  that  genuinely  naodel 
the  physics  of  acoustics  and  convection.  These  components  combine  the  computational  simplicity 
of  FVS  with  the  accuracy  and  stability  of  FDS  and  also  feature  eigenvalues  with  uniform  algebraic 
sign.  This  formulation  eliminates  the  unstable  linear-dependence  problem  in  steady  low-Mach- 
number  flows  and  satisfies  by  design  the  upstream-bias  stability  condition.  As  the  Mach  number 
increases,  the  formulation  smoothly  approaches  and  then  becomes  an  upstream-bias  approximation 
of  the  entire  flux  divergence,  along  one  single  direction. 

4.1  Convection  and  Pressure-Gradient  Components 

The  flux  divergence  can  be  decomposed  into  convection  and  pressure-gradient  components  as 

^  ^  ^  (41) 

dx  dx  dx 

where  P  and  P  respectively  denote  the  convection  and  pressure  fluxes,  defined  as 


f  \ 

'  0  ' 

m 

p 

^  p 

m 

^  =  —  { 
P 

m 

P 

-{E  +  p) 

E  +  p 

0 

p 
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For  supersonic  flows,  the  Euler  eigenvalues  (13),  associated  with  all  have  the  same  algebraic 
sign  and  the  entire  flux  divergence  can  be  upstream  approximated  along  one  single  direction.  For 
subsonic  flows  these  eigenvalues  have  mixed  algebraic  sign  and  an  upstream  approximation  for  the 
flux  divergence  along  one  single  direction  remains  inconsistent  with  the  two-way  propagation  of 
acoustic  waves.  Without  the  pressure  gradient  in  the  momentum  equation,  however,  the  corre¬ 
sponding  flux-jacobian  eigenvalues  all  have  the  same  algebraic  sign^  and  the  resulting  convection 
flux  divergence  can  then  be  upstream  approximated  along  one  single  direction.  The  flux  divergence 
can  thus  be  decomposed  as  the  linear  combination 


dx 


dx  ^  dx 


+ 


dx  . 


0</3<l 


(43) 


where  the  positive  pressure-gradient  partition  function  /?  can  be  chosen  in  such  a  way  that  all  the 
eigenvalues  of  each  of  the  two  components  between  brackets  in  (43)  keep  the  same  algebraic  sign  for 
all  Mach  numbers.  In  this  manner,  these  entire  components  can  be  upstream  approximated  along 
single  directions.  This  choice  for  /3  is  possible  because  the  eigenvalues  of  a  matrix  are  continues 
functions  of  the  matrix  entries^®  and  hence  all  the  eigenvalues  for  the  components  in  (43)  will 
continuously  depend  upon  (3.  The  function  P  will  gradually  increase  toward  1  for  increasing  Mach 
number,  so  that  an  upstream  approximation  for  the  components  in  (43)  smoothly  approaches  and 
then  becomes  an  upstream  approximation  for  the  entire  along  one  single  direction.  Decom¬ 
position  (43)  is  thus  used  for  an  upstream  approximation  of  the  flux  divergence  for  subsonic  and 
supersonic  flows. 

For  low  and  vanishing  Mach  numbers,  decomposition  (43),  however,  is  insufficient  for  an  accu¬ 
rate  upstream  modeling  of  acoustic  waves.  For  a  Mach  number  that  approaches  zero,  the  Euler 
eigenvalues  (13)  can  all  keep  the  same  algebraic  sign  only  if  the  sound  speed  contribution  vanishes, 
which  corresponds  to  a  vanishing  pressure  gradient  contribution  and  hence  P  approaching  zero.® 
But  for  P  approaching  zero,  the  eigenvalues  associated  with  the  components  in  (43)  approach  the 
eigenvalues  of  the  jacobians 


g/®(g) 

dq 


m 


1  m  ,  m 


and 


dfP 

dq 


f  0 

Pp 

0 


1 

2m 


E  +  p  m 
P  P 


0  \ 


Pb 

0 


0 

0 


5  ~  (1  T Pe) 


(44) 


(45) 


Using  the  pressure-derivative  identity  (5)  the  eigenvalues  of  these  jacobians  respectively  are 


^i=^(i+p*)  («) 

=  (47) 

which  certainly  all  keep  the  same  algebraic  sign,  but  for  vanishing  Mach  number  remain  far  less  than 
the  dominant  speed  of  sound  c.  For  low  Mach  numbers,  therefore,  an  upstream  approximation  for 
the  components  in  (43)  would  inaccurately  model  the  physics  of  acoustics.  This  difficulty  is  resolved 
by  further  decomposing  the  pressure  gradient  in  (43)  in  terms  of  a  genuine  acoustic  component, 
for  accurate  upstream  modeling  of  acoustic  waves. 


and 


a;,,  =  o  , 
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4.2  Acoustic  Components 

Following  the  acoustic  equations  (15),  the  flux  divergence  g  can  be  alternatively  decomposed  for 
arbitrary  Mach  numbers  and  corresponding  dependent  variables  p,  m,  and  E  as 

(48) 

dx  dx  dx  dx  ^  dx 

In  this  decomposition,  the  matrices  and  are  defined  as 

/  0  ,  1  ,  0  \  /  0 

Pol  ®  1  Pe 


V 


0 


0 

-P, 

Pe 


0 


) 


\ 


-1 

Pm 

-Pp 

Pe 


0  \ 
0 


0 


(49) 


) 


Heed,  in  particular,  that  no  flux  component  of  f{q)  exists,  of  which  the  Jacobian  equals  The 
eigenUlues  of  the  matrix  have  been  determined  in  closed  form  as 


A? 2  =  0  >  =  -cMpeuI\u\ 


(50) 


which  become  infinitesimal  for  vanishing  M.  The  matrix  can  be  termed  a  “non-linear  coupling” 
matrix  for  it  completes  the  non-linear  coupling  between  convection  and  acoustics  within  (48)  so 
that  the  two  Euler  eigenvalues  Af  3  in  (13)  do  correspond  to  the  sum  of  convection  and  acoustic 
speeds.  Since  decomposition  (48) ’will  be  used  in  the  upstream-bias  formulation  for  small  Mach 
numbers  only  and  considering  that  the  eigenvalues  in  (50)  vanish  for  these  Mach  numbers,  no  need 
exists  to  involve  in  the  upstream-bias  approximation  of  the  flux  Jacobian  (11). 

The  eigenvalues  of  A“  are  exactly  determined  in  closed  form  as 

A?  =  0  ,  A5.3  =  ±c  (51) 


The  matrix  A“,  therefore,  can  be  termed  the  “acoustics”  matrix,  for  its  eigenvalues,  unlike  (46)- 
(47),  equal  the  speed  of  sound  c  for  any  Mach  number.  Despite  its  zero  eigenvalue,  A“  features  a 
complete  set  of  eigenvectors  and  thus  possesses  the  similarity  form 

A“  =  XA“A-^  =XA“+X-^+XA“-A:-^  ,  A“  =  A“+-I-A“-  (52) 


where  A“+  and  A““  respectively  contain  non-negative  and  non-positive  eigenvalues.  The  Euler  flux 
divergence  decomposition  (48)  thus  becomes 

—  =  A:A“+X“^  +  -1-  ^  +  A""^  (53) 

dq 


Since  the  two  acoustics  matrices  at  the  rhs  of  this  expression  respectively  possess  non-  nega¬ 
tive  and  non-positive  eigenvalues,  a  characteristics-bias  approximation  of  these  matrices  involves 
an  upstream  approximation  of  the  first  matrix  and  a  downstream  approximation  of  the  second 
matrix.  These  approximations  naturally  lead  to  the  following  absolute  acoustics-matrix  upstream 

expression 


p,,/c  >  0  > 

0  ,  c  , 

-Pp)pJ{cPe)  .  0  ’ 


Pe^ 

0 


(c^  -Pp)  /c 


dx 


(54) 


Upon  imposing  the  condition  that  the  continuity  and  energy  components  in  this  matrix  product 
should  also  satisfy  the  acoustic-field  results  (18),  leads  to  the  beautifully  simple  result 


=  =  c/^  ,  /  =  identity  matrix  (55) 

'  '  dx  dx 
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which  indicates  for  this  matrix  product  the  equivalence  of  replacing  |A“|  with  the  matrix  cl,  of 
which  all  eigenvalues  approach  +c.  For  acoustic  flows  and  related  dependent  variables  E  and  p,  (55) 
is  exact.  For  non-acoustic  flows  and  related  arbitrary-Mach-number  dependent  variables  E  and  p 
(55)  is  approximate.  This  computationally  advantageous  approximation  on  the  acoustics  upstream- 
bias,  not  on  the  flux  divergence  itself,  therefore,  will  be  used  essentially  in  the  low  subsonic-flow 
Mach  number  regime. 


4.3  Acoustics  Convection  Characteristics-Flux  Divergence 

The  previous  sections  have  shown  that  the  flux  Jacobian  (11)  can  be  equivalently  expressed  as 


dq 


dq  dq\ 


+ 


(1-/3) 


dfP] 


dq  J 


a/9 


(56) 


-t-  AA“" ^  -b  A 
dq 


where  the  first  expression  is  convenient  for  a  characteristics-bias  approximation  for  high-subsonic 
and  supersonic  Mach  numbers  and  the  second  expression  is  convenient  for  low-subsonic  Mach 
numbers. 

A  flux  Jacobian  decomposition  for  all  Mach  numbers  can  thus  be  cast  as  the  linear  combination 


dq 


(l-a)  I 


^  dP 

dq  ^  dq  \ 


(1  -  /3)^  |aA“+A-^  +  XA“-A-i  +  ^ 


with  0  <  a  <  1,  which  leads  to  the  following  acoustics-convection  decomposition  of  the  flux  Jacobian 


^  =  aAA^+X-^  +  aXA^-X-'^  + 
dq 


dp 

[dq 


+  (1  -  ol)P 


dp- 

dq  . 


dq 


-l-aA”" 


(58) 


As  mentioned  in  Section  4.2,  an  upstream  approximation  to  the  Euler-fiux  jacobian  will  be 
developed  by  establishing  upstream  approximations  only  for  the  first  four  terms  in  (58),  where  the 
jacobian  matrix  +  (1  —  is  counted  as  one  term.  The  reason  for  this  coupling  is  that, 

with  reference  to  (43),  the  eigenvalues  of  this  jacobian  matrix  will  all  keep  the  same  algebraic  sign, 
because  (1  —  a)/3  <  /3. 

One  justification  for  this  selective  upstream  formulation  rests  on  the  physical  acoustics  and 
convection  significance  of  these  terms.  For  any  magnitude  of  both  pressure  and  pressure  gradient, 
the  convection  field  uniformly  carries  information  along  streamlines;  hence,  the  entire  can 
receive  an  upstream  bias  along  one  single  direction.  The  matrices  and  XA^'^X  ^ 

account  for  the  bi-modal  propagation  of  acoustic  waves;  these  matrices  axe  thus  used  for  an  acoustics 
upstream  approximation,  for  low  Mach  numbers.  The  pressure  fiux  too  accounts  for  the  bi-modal 
propagation  of  acoustic  waves,  but  in  conjunction  with  As  the  Mach  number  increases  from 

zero,  a  larger  and  larger  fraction  (1  —  cx)l3^^  of  the  pressure  flux  jacobian  can  thus  be  upstreamed 
in  the  same  direction  as  and  along  with  while  (1  ~  a)(l  —  /3)  ^  is  upstreamed  in  the  opposite 
direction.  As  the  Mach  number  increases,  therefore,  a  smaller  and  smaller  fraction  a{XA°‘'^X'^^  4- 
XA^'^X’^^)  of  (XA“+X“^  +XA°'~X~^)  is  upstreamed.  The  upstream-bias  function  a  will  decrease 
and  P  will  increase  as  the  Mach  number  increases,  so  as  to  ensure  physical  significance  of  the  overall 
upstream-bias  approximation  to  the  first  four  terms  in  (58).  The  function  /3,  in  turn,  depends  on 
another  function  5  that  leads  to  simpler  expressions. 


Given  the  algebraic  sign  of  the  eigenvalue  set  of  each  matrix  term  in  (58),  the  associated 
direction  cosines  ae  for  the  upstream-bias  expression  (28)  are 

oi  =  H-1  ,  02  =  -1  ,  03  =  s  =  sgn(tt)  ,  04  =  -s  =  -sgn(u)  ,  05  =  0  (59) 

where  s  =  sgn(o)  denotes  the  algebraic  sign  of  u.  With  (58),  (59),  approximation  (55),  and 
j  _  q;)(2/?  —  1),  the  general  expression  (28)  leads  to  the  acoustics-convection  characteristics 

flux  divergence 


1 


dx  dx 


if 


ei  d  \  (  dq  df’  M” 

dx  dx  \  dx  dx  dx 


where  I  denotes  the  identity  matrix  of  appropriate  size.  In  particular,  the  coupling  of  an  upstream 
approximation  for  {l-a)(3^,  via 03,  with  a  downstream  approximation  for  (1  -a)(l -/3)-^,  via 
04  results  in  an  overall  upstream  approximation  of  the  pressure  gradient,  but  with  variable  weight 
S.  The  operation  count  for  expression  (60)  is  then  comparable  to  that  of  an  FVS  formulation. 
The  terms  in  this  expression,  furthermore,  directly  correspond  to  the  physics  of  acoustics  and 
convection.  For  low  Mach  numbers,  ^  =  0  and  (60)  reduces  to 


df^  _df  d 
dx  dx  dx 


,  f  dq  dP\  1 

j  J 


which  essentially  induces  only  an  acoustics  upstream.  Heed  that  the  components  within  remain 
linearly  independent  of  one  another,  which  avoids  the  linear-dependence  instability  in  the  steady 
low-Mach-number  Euler  equations.  For  supersonic  flow,  a  =  0  and  (5  =  1.  Expression  (60)  in  this 
case  becomes  ^  ^  ^  ,  s  u 

=  sv-fs^l  (62) 

Bx  dx  ax  [  ^  Vdx)  J 

which  corresponds  to  an  upstream  approximation  of  the  entire  Euler  flux  divergence. 


5  Upstream-Bias  Eigenvalues  and  Functions 

The  acoustics-convection  upstream  functions  a  and  S  depend  on  the  Mach  number.  They  are 
determined  by  enforcing  the  upstream  stability  condition  on  the  upstream  matrix  for  (60).  The 
divergence  of  the  characteristics  flux  in  (60)  becomes 


^  =  +  +  (63) 

dx  dx  L  ^  V  dq  j  dx  \ 

The  terms  between  parentheses  collectively  constitute  the  upstream-bias  dissipation  matrix 

(64) 

dq  dq 

Despite  the  formidable  algebraic  complexity  of  A,  all  of  its  eigenvalues  have  been  analytically 
determined  exactly  in  closed  form.  Dividing  through  by  the  speed  of  sound  c,  the  non-dimensional 
form  of  these  eigenvalues  is 

Ai=a-f-M  ,  X2,3  =  a+ U  + 


In  order  to  ensure  physical  signiflcance  for  the  characteristics-bias  flux  within  (60),  hence  for  the 
upstream-bias  approximation  to  decomposition  (58),  the  upstream  bias  functions  a  and  6  will 
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therefore  be  determined  by  forcing  the  upstream-bias  eigenvalues  (65)  to  remain  positive  for  all 
Mach  numbers.  Following  the  considerations  after  (55),  in  particular,  all  these  eigenvalues  must 
converge  to  1  for  vanishing  Mach  number.  Rather  than  prescribing  some  expressions  for  a  and  6 
and  accepting  the  resulting  variations  for  these  eigenvalues,  physically  reasonable  expressions  for 
these  eigenvalues  are  instead  prescribed  and  the  corresponding  functions  for  a  and  6  determined. 

5.1  Eigenvalue  As 

This  eigenvalue  will  correlate  with  the  absolute- value  Euler  eigenvalue  |M  —  1|.  As  a  consequence. 
As  will  vary  between  1  and  1  —  M  foi  0  <  M  <  1  —  cm  and  smoothly  shift  from  1  —  M  to  M  —  1 
within  the  sonic  transition  layer  1  —  eM'^M<l  +  sm-,  where  sm  denotes  a  transition-layer 
parameter;  in  this  work  sm  =  One  expression  for  As  that  remains  smooth  and  meets  these 
requirements  is  the  composite  spline 


HM)  =  { 


1-M 

,  0 

< 

M  <  1  —  £m 

(M-1)2  sm 
Ism  2 

,  1  -  CM 

< 

M  <  1  4-  £m 

(66) 

M-1 

,  1  +  £m 

< 

M 

5.2  Eigenvalue  Ai 

This  eigenvalue  correlates  with  the  non-dimensional  Euler  eigenvalue  M,  but  it  too  has  to  equal  1 
for  M  =  0;  it  then  must  coincide  with  M  for  M  >  1  and  also  remain  greater  than  As,  as  expressed 
through  (66),  for  consistency  with  the  Euler  eigenvalues  (13).  This  condition  in  particular  implies 
Ai  >  ^.  It  thus  follows  that  Ai  will  vary  between  1  and  M  for  0  <  M  <  5  4-eM-  An  expression  for 
Ai  =  Ai(M)  that  remains  smooth  and  meets  all  of  these  requirements  is  the  composite  spline 


1-M 

,  0 

< 

VI 

\  —  £M 

Ai(M)  =  < 

(M-i)2  1+ej^ 

2eM  2 

>  2 

< 

M  < 

^  +  £m 

M 

5  1 

< 

M 

5.3  Upstream-Bias  Functions  a  and  d  and  Eigenvalues 

Prom  Ai  and  As  in  (65),  the  corresponding  expressions  for  both  a  =  a(M)  and  6  —  S{M)  have 
been  exactly  determined  as 

r,(M\  \  (M\  (Ai(M)-As(M))(Ai(M)-A3(M)-1-PbM)  ,  . 

a{M)  =  X,{M)-M  ,  i(M)  = - 1 +p,M(  A,(M)  -  A,(M) ) - 

The  variations  of  the  upstream-bias  functions  a  =  a(M)  and  S  =  S(M)  and  the  corresponding 
eigenvalues  from  (65)  are  presented  in  Figures  2,  3,  respectively. 

Figure  2  indicates  that  the  upstream-bias  functions  as  well  as  their  slopes  remain  continuous  for 
all  Mach  numbers,  with  0<q:<1,  0<J<1  and  a  =  0  for  M  >  £M)  <5(M)  =  1,  for  M  >  \+£m- 

As  5  =  6{M)  rises,  the  upstream-bias  contribution  from  the  acoustics  matrix  decreases  rapidly, 
reducing  to  less  than  25%  of  its  maximum  at  M  =  0.39  with,  therefore,  concurrent  reduction  of  the 
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effect  of  the  acoustic-flow  formula  (55).  The  variation  of  J  =  S(M)  shows  that  the  pressure-gradient 
contribution  to  this  upstream-bias  formulation  increases  monotonically,  while  remaining  less  than 
25%  of  its  maximum,  for  0  <  M  <  0.7.  When  6(M}  =  1  for  supersonic  Mach  numbers,  the 
entire  pressure-gradient  is  upstreamed  with  the  same  weight  as  in  the  convection  flux,  in  complete 
agreement  with  the  physical  mono-axial  wave  propagation  within  supersonic  flows. 


Figure  3  shows  that  within  0  <  M  <  1  -1-  £m,  the  eigenvalues  Ai,  A2,  A3  smoothly  approach  1 
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for  vanishing  M,  indicating  a  physically  consistent  upstream-bias  approximation  of  the  acoustic 
equations  embedded  within  the  Euler  equations.  For  M  >  I  +  em,  the  eigenvalues  (65)  respectively 
coincide  with  the  Euler  flux  Jacobian  eigenvalues  M,  M  +  1,  M  —  1,  which  corresponds  to  an 
upstream-bias  approximation  of  the  entire  flux  vector,  for  supersonic  flows. 


6  Finite  Element  Weak  Statement 


With  reference  to  (21),  the  divergence  (60)  of  the  characteristics-bias  flux  leads  to  the  following 
characteristics-bias  integral  statement 


dx 


dp  ,dp 

+  s-T — 1- 

ox  ox 


dQ  =  0 


(69) 


An  integration  by  parts  of  the  upstream-bias  expression  then  generates  the  weak  statement 


dn  =  o 


(70) 


where  the  surface  integral  on  dCl  corresponding  to  the  upstream-bias  expression  vanishes  because 
of  the  boundary  condition  ipixQn)  =  0,  imposed  to  eliminate  unnecessary  bormdary  upstream  bias. 
The  discrete  equations  then  result  from  a  flnite  element  discretization  of  this  weak  statement. 


6.1  Galerkin  Finite  Element  Equations 

The  finite  element  weak  statement^’^’®  associated  with  (70)  is 

dn  =  0  (71) 

where  superscript  “h”  signifies  spatial  discrete  approximation.  The  approximation  exists  on  a 
partition  C  Q,  of  Cl.  This  partition  Cl^  has  its  boundary  nodes  on  the  boundary  dCl  of  O  and 

results  from  the  union  of  Ne  non-overlapping  elements  Cle,Cl^  =  U^i  ^e-  For  mesh  nodes  within 
there  exist  clusters  of  “master”  elements  each  comprising  only  those  adjacent  elements 
that  share  a  mesh  node  Xi,  which  implies  existence  of  exactly  N  master  elements.  As  Figure  4 
shows,  on  each  master  element  the  discrete  test  function  =  Wi  =  Wi  (x),  1  <  i  <  N,  will 
coincide  with  the  “pyramid”  basis  function  with  compact  support  on  Clf'.  Such  a  function  equals 
one  at  node  Xi,  zero  at  all  other  mesh  nodes,  and  also  identically  vanishes  both  on  the  boundary 
segments  of  fip  not  containing  Sj,  and  elsewhere  within  the  computational  domain  outside  Clf'. 

Note  that  represents  a  “finite  volume”  as  used  in  finite  volume  schemes,  which,  however, 
do  not  employ  pyramid  test  functions.  The  following  developments  are  based  on  a  linear  pyramid 
test  function  Wi,  which  can  be  expressed  as 


Wi{x)  =  < 


X  -  Xj-l 

Ax,  1 

®  2 

Xj-i-l  -  X 


(72) 


The  discrete  solution  at  each  time  t  assumes  the  form  of  the  following  linear  combination 

N 

(x,  t)  =  Y^  Wj  (x)  •  q^  {xj,t) 
i=i 


(73) 


Figure  4:  Master  Element  fif  and  Test  Function  Wi  =  Wi{x) 

of  nodal  solution  values  and  trial  functions,  which  coincide  with  the  test  functions  Wj{x)  for 
a  Galerkin  formulation.  Similarly,  the  source  </>  =  4>(x,q{x,t))  and  fluxes  /  =  f{q{x,t)), 

jq  ^  /9(q(a;,t))  and /P  =  /^(^(x,*))  are  discretized  through  the  group  expressions 

<f)^  {x,  t)  =  ^  Wj  (x)  •  (j)  (xj,  q^  {Xj,  t))  ,  {x,  t)  =  52  ‘^3  ^ 


/«'*  (x,  t)  =  ^  Wj  (x)  •  /«  (q^  (x,-,  t))  ,  /P'  (x,  t)  =  j2  (^)  •  Z"’  *)) 

jf=i 

The  notation  for  the  discrete  nodal  variable  and  fluxes  is  then  simplified  as  qj{t)  =  q^{xj,t), 
fj  =  f]  =  /«'*  {xj,t),  /J  =  /p'  {Xj,t)  and  expansions  (73)-(74)  are 

then  inserted  into  (71),  which  yields  the  discrete  finite  element  weak  statement 

for  1  <  «  <  iV,  with  £  set  equal  to  a  reference  length  within  each  element,  typically  a  measure 
of  the  element  size.  While  an  expression  like  (73)  for  ip^,  a^,  ,  s^,  and  can  be  directly 
accommodated  within  (75),  each  of  these  variables  in  this  study  has  been  set  equal  to  a  piece 
wise  constant  for  computational  simplicity,  one  centroidal  constant  value  per  element.  Since  the 
test  and  trial  functions  Wi  are  prescribed  functions  of  x,  the  spatial  integrations  in  (75)  are  then 
exactly  carried  out,  which  transforms  (75)  into  a  system  of  ordinary  differential  equations  (ODE) 
in  continuum  time  for  determining  at  each  time  level  t  the  unknown  nodal  values  q^  {xj,t). 

For  a  clear  comparison  between  traditional  finite  difference/volume  schemes®  and  the  acoustics- 
convection  finite  element  algorithm  (75),  at  any  interior  node  “i”  of  the  representative  grid  in 
Figure  4,  equations  (75)  with  e'^  =  (Axi+i)/2  can  be  equivalently  recast  in  difference  notation  as 


=  -{xpac).!  [qi  -  qi-i)  +  (ft+i  -  qi)  + 


( {fi  -  fi-i)  •  (i  +  («V’)i_i)  +  {fi+i  -  fi)  •  (i  )  + 

-I  ( (/f  -  /f_x)  •  (i  +  (^#)i-|)  +  (/IVi  -  /f)  •  (i  -  ) 


(76) 
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which  uniquely  couples  several  time  derivatives  at  each  node  “z”  and  features  a  linear  combination 
of  two-point  upstream  and  downstream  flux  differences.  In  these  finite  element  equations,  the 
values  of  the  controller  -0^  determines  the  combination  weights  of  the  downstream  and  upstream 
expressions,  and  since  0^  remains  non-negative,  these  equations  induce  the  appropriate  upstream 
bias,  since  the  upstream  weight  1  +  always  exceeds  the  downstream  weight  1  —  As  a 

result,  the  finite  element  weak  statement  (75)  generates  consistent  vaxiable-upstream-bias  discrete 
equations  that  correspond  to  an  upstream-bias  discretizations  for  the  original  Euler  system  (1), 
within  a  compact  block  tri-diagonal  matrix  statement. 

For  smooth  solutions,  these  equations  will  still  couple  upstream  and  downstream  points  even 
for  supersonic  flows.  The  potential  objection  that  one  such  algorithm  would  violate  the  physics 
of  mono-directional  wave  propagation  for  supersonic  flows  is  easily  addressed  with  Courant’s  and 
Hilbert’s  classical  developments^^  for  non-linear  hyperbolic  systems.  They  in  fact  concluded  that 
while  waves  propagate  along  characteristics,  smooth  solutions  can  be  expanded  in  Taylor’s  series 
within  arbitrary  regions  encircling  any  given  point  and  along  any  direction  radiating  upstream  or 
downstream  from  the  point. 

For  a  closer  comparison  with  upwind  finite- volume  schemes^,  the  finite  element  equations  (76) 
can  be  rearranged  to  generate  the  “numerical  flux” 


_  fi  +  fi+l 


-^i+h 


J  -  „)  +  5±i  -  ff)  +  rp.  -  z?) 


(77) 


which  corresponds  to  the  discrete  counterpart  of  the  characteristics-bias  flux  within  (60) .  By  virtue 
of  this  numerical  flux,  equations  (76)  are  recast  as 


(dqj-i  dgj 

6  \  dt  dt 


4^i~l 


+ 


dqi+1 

6  \dt  dt 


(78) 


which  shows  that  the  finite  element  weak  statement  (75)  naturally  leads  to  a  discretely  conservative 
algorithm. 


6.2  Boundary  Equations  and  Pressure  Boundary  Condition 

The  integral  statement  (70)  directly  yields  a  set  of  consistent  boundary  differential  equations,  for 
both  imconstrained  boundary  variables  and  for  pressure,  to  enforce  a  pressmre  boundary  condition 
at  a  subsonic  outlet.  These  equations  do  not  require  any  algebraic  extrapolation  of  variables,  but 
rather  couple  the  time  derivatives  of  boundary-  and  interior-node  variables  within  the  boundary 
cell. 

For  the  linear  elements  in  this  study,  let  N  and  IV  —  1  denote  the  nodes  within  the  oulet  bound¬ 
ary  element,  with  N  corresponding  to  the  outlet  node.  For  the  discrete  finite  element  equation 
associated  with  boundary  node  Xjv,  the  controller  ^  and  test  function  w  satisfy  the  conditions 
^  =  0  and  w{xn-i)  =  0,  w{xn)  =  1. 

The  boundary  differential  equation  from  (75)  corresponding  to  an  outlet  node  becomes 

This  equation  directly  couples  the  time  derivatives  of  the  solution  q  at  the  adjacent  boundary  and 
interior  nodes  Xpf  and  Xjv-i-  A  similar  equation  is  then  obtained  at  an  inlet,  mutatis  mutandis. 
Furthermore,  no  upstream  bias  is  necessary  within  a  boundary  equation,  hence  0  =  0,  because  as 
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(79)  shows,  this  finite  element  boundary  equation  directly  yields  an  upwind  approximation  for  the 
divergence  of  /. 

Concerning  the  pressure  outlet  boundary  condition,  this  is  natmrally  enforced  within  the  surface 
integral  that  emerges  in  the  momentum-equation  weak  statement.  The  convection  and  pressure 
flux  decomposition 

m  =  nq)+F{q)  (80) 

is  first  inserted  into  the  non-discrete  integral  statement  (70);  subsequent  integration  by  parts  of 
the  pressure  gradient  therein,  generates  the  weak  statement 


w{x)f^ 


■]  X=XJV 


=  0 


X—XN-\ 


(81) 


A  subsequent  linear  finite  element  discretization  of  (81)  yields 

=  -  (^«  -  -  (  PA.  (82) 

In  this  equation,  denotes  the  outlet-node  pressure,  as  calculated  through  the  equation  of  state 
(8),  whereas,  quite  significantly,  can  correspond  to  the  specified  outlet  pressure  boundary 
condition.  This  strategy  for  imposing  an  outlet  pressure  boundary  condition  remains  intrinsically 
stable.  Suppose,  for  instance,  that  some  numerical  perturbation  forces  to  decrease  below  the 
imposed  f^ut-  I’l  case  the  outlet  boundary  equation  (82)  induces  a  negative  time  rate  of 
change  for  which  leads  to  a  corresponding  reduction  in  rrij^.  From  the  equation  of  state  (8), 
this  reduction  then  leads  to  an  increase  in  /^,  which  corresponds  to  a  stable  restoration  of  the 
imposed  pressure  condition.  A  similar  conclusion  on  the  stability  of  (82)  is  achieved  by  considering 
a  perturbation  increase  in  /^.  The  results  in  Section  8  confirm  tha  accuracy  and  stability  of  this 
pressure  boundary  condition  procedure. 


6.3  Discrete  Upstream-Bias  Controller 

This  section  introduces  a  new  upstream-bias  controller  This  controller  varies  in  the  range 
0  <  V’min  <  <  ^max  <  1  and  controls  within  (76)  the  amount  of  induced  upstream-bias,  hence 

artificial  difllusion.  By  design,  =  0  corresponds  to  a  classical  centered  discretization,  hence  no 
induced  diflFusion;  =  1,  instead,  corresponds  to  a  fully  upwind  formulation,  hence  maximum 
diffusion. 

Denote  then  with  ipi  the  numerical  value  of  the  controller  at  the  representative  node  “z”.  By 
analogy  with  (73),  the  discrete  controller  is  cast  as  the  following  linear  expansion 

N  N 

{x,  t)  =  Yl  >  *)  =  (83) 

j=l  .7=1 


and  the  centroidal  evaluation  of  this  expression  within  each  element  then  yields 

,  ipi  +  ^i+i 

M  =  — 2 — 


(84) 


In  regions  of  smooth  flow  approaches  V’min  for  a  local  reduction  of  upstream-bias  diflfusion;  in 
region  of  discontinuous  solution  approaches  ipmax^  for  an  essentially  non-oscillatory  resolution  of 
local  discontinuities.  The  controller  will  thus  correlate  with  a  local  measure  (pi  of  slope  discontinuity. 
For  a  smooth  variation,  ■ipi  is  expressed  as  a  function  of  the  measure  (pi  as  the  composite  spline 
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ipi 


t 


< 


i’max  + 


V’min 

V’max  ^min 
{‘Pd  -  V>c? 


-  {Pd  -  ^Pc)  <^1  +  3  {(pr>  +  Pc)  p\  -  2</5f 


Pi  <  Pc 


Pc  <  Pi  <  Pd 


1^  V^max  j 

(85) 

where  0  <  <  1  and  the  positive  and  respectively  denote  threshold  continuity  and 

discontinuity  measures.  Figure  5  shows  the  smooth  variation  of  this  type  of  spline  controller. 


Figure  5:  Variation  of  Controller  ipi 

The  implementation  of  thus  requires  a  set  of  points  within  where  the  slopes  of  the 
approximate  solution  are  generally  discontinuous.  For  the  finite  element  approximation  (73),  this 
set  of  points  is  the  set  of  finite-element  side  nodes  shared  by  distinct  finite  elements  in  for 
the  continuous  expansion  (73)  changes  approximating  polynomial  from  element  to  element,  which 
implies  that  solution  slopes,  are  generally  discontinuous  at  element  side  nodes.  This  type  of  slope 
discontinuity  is  depicted  in  Figure  6  via  the  local  normal  unit  vectors  and  respectively  to 
the  left  and  right  of  the  slope-discontinuity  point  P,  where  there  exist  two  distinct  normal  unit 
vectors,  one  for  each  of  two  elements  sharing  the  node. 

With  reference  to  Figure  6,  the  magnitude  ||n^  — n^||  of  the  vector  difference  n^  —  n^  becomes 
proportional  to  a  bounded  measure  of  local  slope  discontinuity.  If  the  graph  slope  is  continuous 
at  P,  then  coincides  with  and  ||n^  —  n^||  vanishes.  On  the  other  hand,  when  a  slope 
discontinuity  exists  at  P,  as  shown  in  the  figure,  ||n^  —  n^||  varies  between  0  and  2  depending 
on  the  magnitude  of  the  slope  jump.  A  positive  measure  of  slope  discontinuity  that  vanishes  for 
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Figure  6:  Slope  Discontinuities  and  Local  Unit  Vectors 


continuous  slopes,  remains  bounded  and  strictly  varies  between  0  and  1  can  thus  be  defined  as 

-  n^\\x=xi  (86) 

By  virtue  of  the  law  of  cosines,  the  local  measure  ipi  also  equals 


where  9  denotes  the  angle  between  the  unit  vectors  in  Figure  6. 

With  reference  to  Figure  7  and  (87),  specific  numerical  values  for  (pp,  V’min  and  V’max 
can  be  easily  established.  At  a  point  of  solution  smoothness,  like  point  i  —  1  in  the  figure, 
will  be  parallel  to  n^,  hence  9  =  0°  which  fi'om  (87)  leads  to  =0.  At  a  shock,  instead,  9 
can  become  greater  than  90°,  as  shown  in  the  figure  for  point  i.  The  threshold  9  =  90°  is  thus 
selected  for  ip^,  which  from  (87)  leads  to  =  l/\/2.  Numerous  numerical  experiments  with  the 
acoustics-convection  algorithm  have  indicated  that  a  minimal  amount  of  “background”  upstream 
bias  is  necessary  for  convergence;  this  finding  is  not  surprising,  since  the  formulation  is  essentially 
centered,  hence  devoid  of  any  upstream-bias  dissipation  for  =  0.  Hence,  tpmm  >  0,  with  typical 
numerical  values  in  the  range  |  <  tpmin  <  5-  Concerning  V’maxi  a  relation  with  readily  follows 
from  the  requirement  that  in  the  neighborhood  of  a  shock  the  maximum  upstream  bias  can  at 
most  correspond  to  a  fully  upwind  algorithm,  for  an  essentially  non-oscillatory  capturing  of  shock 
waves.  Hence,  from  (84),  '0t+i  ^  1  with  i  -H  5  denoting  the  centroid  of  a  finite  element  (  cell  ) 
that  supports  a  shock.  For  a  typical  case  of  a  shock  captured  within  at  least  two  cells,  as  shown 
in  Figure  7,  (85)  leads  to  ipi  =  V’max  and  ipt+i  V’min-  From  (84),  therefore 


'Pmax  "I"  '0min 
2 


<  1  'Pmax  ^  2  Ipmm 


(88) 


which  linearly  decreases  as  a  function  of  increasing  ipmin-  The  specific  objective  of  letting  vary 
as  the  solution  evolves  is  to  minimize  induced  upstream-bias  dissipation  for  maximum  accuracy 
within  the  prescribed  computational  stencil.  As  its  distinguishing  design  feature,  the  acoustics- 
convection  upstream  resolution  algorithm  nevertheless  remains  an  authentic  characteristics-bias 
formulation  for  any  ip^  with  ipmin  ^  ^  V^max- 


X 


Figure  7:  Local  Unit  Vectors  at  a  Shock 

The  general  expression  of  (pi  corresponding  to  a  scalar  component  of  directly  derives 
from  the  finite  element  expansion  (73),  which  can  be  expressed  in  synthetic  implicit  form  as 
=  9c  ~  9c  (®)*)  =  0-  Hence,  a  normal  unit  vector  n  can  be  cast  at  eaoh  time  level 
t  as  n  =  gradF(g^,a;,t)/||gradF||,  where  the  vector  operator  “grad”  encompasses  the  dependent 
variable  q^ .  The  expression  for  the  corresponding  ipi  at  time  level  t  and  point  Xi  then  becomes 


where  the  partial  derivatives  are  determined  through  the  finite  element  expansion  (73),  and  where 
superscripts  L  and  R  indicate  evaluation  within  the  elements  respectively  to  the  left  and  right  of 
X  =  Xi.  The  form  of  (89)  at  node  i  of  a  uniform  grid  is 


where  the  denominator  never  vanishes.  This  expression,  furthermore,  remains  bounded  and  dif¬ 
ferentiable  for  arbitrary  nodal  values  of  q^.  For  the  sole  purpose  of  determining  the  order  of  this 
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expression  with  respect  to  Ax,  for  a  smooth  solution  over  any  two  contiguous  elements,  the  discrete 
solution  values  “  1  <  J  <  *  +  !>  over  these  elements,  can  be  considered  as  the  nodal  values  of  a 
single  auxiliary  continuous  functions  qc{x,t):  a  Lagrangian,  trigonometric,  or  other  interpolant  of 
the  Qcj ’s  over  both  elements.  With  this  consideration,  the  Taylor’s  series  expansion  of  (90)  yields 

1  +  {q'cixi,  ^  ^  ^  (91) 

where  superscript  prime  indicates  differentiation  with  respect  to  x  and  /C  denotes  the  local  curva¬ 
ture.  This  expansion  reveals  that  tpi  decreases  for  vanishing  Ax.  Even  for  large  slopes,  furthermore, 
tpi  remains  of  order  Ax  in  regions  of  small  curvature.  Only  when  both  curvature  and  slope  drasti¬ 
cally  rise,  e.g.  at  a  shock,  will  ipi  increase,  which  precisely  corresponds  to  the  desired  behavior. 


= 


Ax 

~77Z~^~2 


+ 


O  (Ax2)  =  K 


7  Implicit  Runge  Kutta  Time  Integration 

The  finite  element  equations  (76)  along  with  appropriate  boundary  equations  and  conditions,  de¬ 
lineated  in  Sections  6.2  and  8,  can  be  abridged  as  the  non-linear  ODE  system 

M^^  =  F{t,Q{t))  (92) 

where  corresponds  to  the  coupling  of  time  derivatives  in  (76),  and  F  {t,  Q{t))  represents  the 

remaining  terms  in  (76).  The  numerical  time  integration  of  (92)  in  this  study  takes  place  through 
a  new  class  of  two-stage  diagonally  implicit  Runge-Kutta  algorithms^  (IRK2)  expressed  as 

Qn+l  —  Qn  —  biKi+  b2K2 

A4Ki  =  At  ■  F  {tn  +  ciAt,Qn  +  aiiKi) 

MK2  —  At  ■  F  {tfi  -l-  C2  At,  Qn  +  ci2iKi  -l-  (I22K2)  (93) 

where  n  now  denotes  a  discrete  time  station  and  61,  h2,  ci,  C2,  an,  021,  and  022  indicate  constant 
Runge-Kutta  coefficients,  subject  to  the  constraints  Cj  =  X)i=i  =  1-  The  coefficients 

for  second  order  accuracy  are  listed  in  the  Table  1.  With  these  coefficients,  in  particular,  algorithm 

Table  1:  Runge-Kutta  Coefficients 


bi  62  an  “21  022 

IRK2 

3-v^  l-H\/3  3-y/3  ^  v/3-1 

4  4  6  2 

(93)  becomes  absolutely  stable  for  arbitrary  stiff  non-linear  dissipative  ODE  systems.^’^^ 

This  algorithm  is  implicit  because  the  entries  in  the  arrays  Ki  and  K2  remain  coupled  and  are 
then  computed  by  solving  algebraic  systems.  Diagonally  implicit  signifies  that  Ki  is  determined 
independently  of  ifa-  Thus,  given  the  solution  Qn  at  time  Ki  is  computed  first,  followed  by 
K2.  The  solution  Qn+i  is  then  determined  by  way  of  the  first  expression  in  (93). 

The  terminal  numerical  solution  is  then  determined  using  Newton’s  method,  which  for  the 
implicit  fully-coupled  computation  of  the  IRK2  arrays  K,,  1  <  f  <  2,  is  cast  as 


M  -  a^At  -  JCf)  =  At  F  (t„  -F  a  At,  Qf)  -  MK[ 
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Qf^Qn  +  anKf-hai2KP  (94) 

where  Oy  =  0  for  j  >  i,  p  is  the  iteration  index,  and  Kf  —  Ki  for  i  =  2\  for  linear  finite  elements, 
the  Jacobian 

/dF\^ 

JiiQ)  =  M-  auM  (95) 

then  becomes  a  block  tri-diagonal  matrix.  For  all  the  results  documented  in  the  discussion  section, 
the  initial  estimate  is  set  equal  to  the  zero  array,  while  only  one  iteration  is  executed  for  (94) 
within  each  time  interval.  In  this  mode,  Newton’s  iteration  becomes  akin  to  a  classical  direct 
linearized  implicit  solver. 


8  Computational  Results 


The  computational  results  have  validated  the  accvuracy  and  essential-monotonicity  performance 
of  the  acoustics-convection  upstream  resolution  algorithm  for  transient  and  steady  smooth  and 
shocked  mixed  subsonic  /  supersonic  flows.  The  algorithm  has  generated  essentially  non-oscillatory 
results  that  automatically  preserve  a  constant  total  enthalpy  as  well  as  smoothness  of  both  enthalpy 
and  linear  momentum  across  steady  normal  shocks.  These  results  reflect  available  exact  solutions 
and  numerical  results  independently  generated^^  using  van  Leer’s  and  Roe’s  schemes.  The  bench¬ 
marks  in  this  section  cover  a  total  of  5  different  perfect-gas  flows  encompassing  flows  within:  a  shock 
tube,  a  convergent-divergent  (  DeLaval )  nozzle  and  a  steeply  diverging  nozzle.  The  corresponding 
spatial  computational  domain  fl  for  all  the  results  presented  is  defined  as:  fi  =  [a,  b]  =  [0, 1] , 
uniformly  discretized  into  100  linear  finite  elements,  hence  Ax  =  0.01.  For  each  benchmark,  the 
calculations  proceeded  with  a  prescribed  constant  maximum  Courant  number  Cmax  defined  as 


Cmax  =  max{|tt  -1-  c(,  \u  -  c|,  c} 


At 

Ax 


(96) 


Given  Ax  and  Cmax  for  each  benchmark,  the  corresponding  At  was  thus  determined  as 


At  = 


_ Cmax^^ _ 

max{|u  4-  c|,  |it  —  c|,  c} 


(97) 


As  detailed  in  Section  6.3,  the  upstream-bias  controller  uses  one  scalar  component  of  the  dependent 
variable  q.  In  this  study,  the  algorithm  has  employed  total  energy  E  to  calculate  0. 

All  the  solutions  in  these  validations  are  presented  in  non  dimensional  form,  with  density  p, 
pressure  p,  energy  E  and  enthalpy  JI  made  dimensionless  through  their  respective  inlet  stagnation 
( total )  values.  The  non-dimensional  speed  is  obtained  by  way  of  the  stagnation  speed  of  sound 

divided  by  ./y.  The  reference  speed  Ur  thus  becomes  Ur  =  y^7Ptoti„/ptoti„/v7  = 

Linear  momentum  is  then  made  dimensionless  through  total  inlet  density  and  reference  speed. 


8.1  Shock- Tube  Flow 

This  benchmark  consists  in  determining  the  gas  flow  that  evolves  from  a  rest  state  within  a  straight 
tube.  The  tube  is  initially  divided  into  two  chambers  separated  by  an  impermeable  diaphragm 
placed  at  the  midpoint  on  the  tube  axis.  The  non-dimensional  initial  conditions  for  the  gas  in  each 
of  the  two  chambers  are 

p=1.00  ,  m  =  0.00  ,  E  =  2.50  ,  0.0  <  x  <  0.5 
p  =  0.125  ,  m  =  0.00  ,  E  =  0.25  ,  0.5  <  x  <  1.0  (98) 
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The  diaphragm  ruptures  at  <  =  0  and  the  solution  corresponding  to  t  =  0.14152  is  sought.  At 
this  time  station,  the  exact  solution  features  a  normal  shock  centered  at  x  =  0.75,  for  each  of  the 
components  of  the  dependent  variables  in  g;  the  distribution  of  density  p  also  develops  a  contact 
discontinuity  centered  at  a:  =  0.62.  Figures  8  a)-c)  present  three  density  distributions  as  obtained 
with  the  acoustics-convection  upstream  resolution  algorithm  for  Cmax  =  l-O,  but  with  a  constant 
upstream-bias  controller  ip. 
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Figure  8:  Density  Upstream  Bias:  a)  100%,  b)  75%,  c)  50%,  d)  25%-50%  Controlled 


The  numerical  values  for  ip  for  these  three  solutions  are  ip  =  1.0,  ip  =  0.75,  ip  =  0.5,  respectively 
corresponding  to  100%,  75%  and  50%  upstream  bias.  The  figures  show  that  a  decrease  in  upstream- 
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bias,  hence  associated  artificial  diffusion,  corresponds  to  an  expected  increase  in  solution  resolution. 
The  solution  in  Figure  1  c),  for  V’  =  0.5,  already  displays  correct  solution  features  with  contact 
discontinuity  and  normal  shock  centered  at  the  exact  locations,  even  though  the  sharp  decrease  in 
density  modehng  the  contact  discontinuity  appears  somewhat  diffused.  Note,  however,  that  this 
solution  remains  essentially  non-oscillatory  without  employing  a  fully-upwind  discretization  with 
MUSCL-type  extrapolation  of  variables. 

Figures  8  d)-9  b)  present  the  solution  generated  with  a  variable  controller  ip,  with  0.25  <ip  < 
0.62  and  Cmax  =  1-0.  This  solution  remains  essentially  non-oscillatory  throughout  the  computa¬ 
tional  domain.  The  contact  discontinuity  in  Figure  8  d)  is  now  resolved  over  about  five  nodes, 
with  increased  overall  solution  sharpness.  The  normal  shock  is  captured  over  two  nodes  and  the 
two  plateaus  juxtaposed  the  contact  discontinuity  remain  essentially  fiat.  Figure  9  a)  presents  the 
corresponding  distribution  of  total  energy  E,  which  reflects  the  features  in  the  density  solution. 


Figure  9:  Controlled  Upstream:  a)  Energy  and  Controller,  b)  Mach  Number,  a  and  6 


This  figure  also  shows  the  associated  variation  of  ip,  indicating  that  ip  remains  close  to  its 
minimum  over  the  smooth  parts  of  the  solution.  Only  at  slope  discontinuities  does  ip  increase, 
following  its  design  features  documented  in  Section  6.3.  Therefore,  ip  increases  marginally,  at  the 
expansion  extremities  and  contact  discontinuity,  and  more  markedly  at  the  normal  shock.  At  these 
locations  the  energy  slopes  abruptly  change.  In  particular,  the  larger  increases  of  ip  remain  localized 
at  the  shock  region,  where  it  is  precisely  needed  for  an  essentially  non-oscillatory  capturing  of  the 
shock.  This  shock  is  captured  with  ip  reaching  0.5,  which  corresponds  to  just  50%  upstream  bias. 
Away  from  the  shock  and  other  slope  discontinuities,  ip  =  0.25  which  corresponds  to  a  mere  25% 
upstream  bias.  This  solution  is  therefore  achieved  with  an  essentially  centered  discretization,  which 
leads  to  the  conclusion  that  a  uniformly  fully  upwind  formulation  is  not  strictly  necessary  within 
a  characteristics-bias  algorithm  to  generate  an  essentially  sharp  and  non-oscillatory  solution. 
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Figure  9  b)  presents  the  distributions  of  Mach  Number  M  and  associated  acoustics  and  con¬ 
vection  upstream-bias  functions  a  and  S.  The  solution  for  M  correlates  with  that  for  p  and  E  with 
essentially  flat  plateaus  and  sharp  normal  shock.  The  distribution  of  a  indicates  that  the  acous¬ 
tics  upstream-bias  induced  via  the  absolute  acoustics  matrix  (55)  remains  significant,  at  a  level 
greater  than  30%,  only  for  M  <  0.3.  Without  this  acoustics  upstream  bias,  essential  monotonicity 
is  lost.  For  these  Mach  numbers,  (5  =  0,  which  corresponds  to  a  centered  approximation  of  the 
pressure  gradient,  for  any  ‘tp.  For  M  >  0.3,  a  decreases  sharply,  accompanied  by  a  corresponding 
rise  in  6.  In  particular,  a  =  0  for  the  Mach  numbers  corresponding  to  the  contact  discontinuity, 
which  indicates  that  the  absolute  acoustics-matrix  upstream-bias  term  acdq/dx  plays  no  local  role 
in  the  calculation  of  this  discontinuity.  With  the  calculated  values  of  the  characteristics-bias 
formulation  for  these  Mach  numbers  is  realized  only  through  the  convection  and  pressure  gradient 
components  intrinsic  to  the  Euler  flux  divergence.  The  aesthetic  appearance  of  the  calculated  con¬ 
tact  discontinuity,  therefore,  is  essentially  due  to  the  linear  flux  approximation  employed,  which 
only  uses  two  nodes  per  cell  without  extrapolation  of  variables.  As  M  increases  towards  its  higher 
plateau,  5  reaches  its  peak  numerical  value  of  about  0.75,  which,  with  xp  =  0.35,  corresponds  to  a 
mere  26.3%  upstream-bias  in  the  approximation  of  the  pressure  gradient.  This  level  of  upstream 
bias  further  decreases  toward  the  shock,  settling  to  a  level  of  15%  with  6  ci  0.3  and  xp  c;  0.5  at  the 
shock. 

Figures  10  a)-b)  display  the  associated  distributions  for  speed  and  static  pressure.  The  plateaus 
in  these  distributions  remain  flat  and  unperturbed  by  the  density  contact  discontinuity,  with  sharp 
shocks  captured  over  about  two  nodes. 
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Figure  10:  Controlled  Upstream:  a)  Speed,  b)  Pressure 
8.2  Flows  in  a  Converging-Diverging  Nozzle 

These  benchmarks  test  the  capability  of  the  algorithm  to  calculate  steady  isentropic  and  shocked 
flows  than  contain  a  low-Mach  number  subsonic  region.  The  nozzle  cross-section  area  distribution 
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remains  continuous  with  continuous  slopes,  but  contains  a  discontinuous  throat  curvature^®,  as 
shown  in  Figure  11.  This  feature  will  induce  a  nozzle-throat  discontinuous  cmvature  in  the  flow 
variables  even  for  an  isentropic  flow,  which  makes  this  benchmark  particularly  useful  to  assess 
algorithm  resolution. 

The  nozzle  area  ratio  distribution  for  these  benchmarks  is 
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a: 

Figure  11:  Variation  of  Area  Ratio  A(x)/A» 

The  initial  conditions  for  the  gas  correspond  to  an  M  =  0.23954  uniform  state,  leading  to  the 
following  initial  numerical  values  for  p,  m,  and  F  throughout  the  nozzle 

p  =  0.97188  ,  m  =  0.27389  ,  F?  =  2.44072  (101) 

The  inlet  remains  subsonic  and  the  physically  admissible  boundary  conditions  specified  at  the  inlet 
are  constant  Dirichlet  conditions  on  density  p  and  total  energy  E,  equal  to  the  initial  conditions. 
An  outlet  boundary  condition  is  also  imposed  at  a  subsonic  outlet,  as  discussed  in  Section  6.2. 

8.2.1  Isentropic  Supersonic  Flow 

For  the  given  initial  conditions,  the  outlet  is  temporarily  subsonic  and  a  pressure  boundary  con¬ 
dition  is  imposed.  This  corresponds  to  the  jump  decrease  in  static  pressiure:  p/ptoti„  =  0.16017, 
pertinent  to  the  terminal  steady-state  outlet  supersonic  Mach  number  M  =  1.85413.  The  algo¬ 
rithm  monitors  the  outlet  Mach  number  during  the  calculations  and  when  it  exceeds  one,  the  outlet 


30 


Joe  lannelli 


pressure  boundary  condition  is  released,  allowing  the  calculations  to  proceed  toward  a  steady  state 
without  any  further  outlet  boundary  condition. 

Figures  12  a)-b)  present  the  convergence  rate  and  steady-state  speed  and  Figures  13  a)-b)  show 
the  distributions  of  density,  linear  momentum,  pressure  and  enthalpy. 


Figure  12:  Isentropic  Flow:  a)  Convergence  Rate,  b)  Speed 


Figure  13:  Isentropic  Flow:  a)  Density  and  Momentum,  b)  Pressure  and  Enthalpy 
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The  acoustics-convection  upstream  resolution  algorithm  generated  a  steady  state  with  a  reduc¬ 
tion  of  the  maximum  residual  by  13  orders  of  magnitude,  down  to  machine  zero  in  about  25  time 
steps  with  Cinax  =  400.  The  calculated  speed  reflects  the  exact  solution,  indicated  with  a  solid 
line.  This  distribution  clearly  shows  the  sonic-point  curvature  discontinuity,  which  remains  devoid 
of  any  unphysical  expansion  shock  and  follows  the  exact  solution. 

The  distributions  of  density,  linear  momentum,  pressure  and  enthalpy  in  Figures  13  a)-b)  also 
reflect  the  corresponding  exact  solutions  with  the  curvature  discontinuity  clearly  resolved.  The  algo¬ 
rithm  has  also  correctly  held  constant  the  computed  enthalpy,  which  satisfies  the  steady-adiabatic- 
flow  constant-enthalpy  condition.  Figure  14  a)  presents  the  distributions  of  total  energy  E,  which 
visually  coincides  with  the  exact  solution,  and  upstream  bias  controller  ip,  with  0.25  <  ip  <  1.0. 
Since  the  solution  does  not  contain  any  shock,  ip  varies  modestly,  with  ip  ~  0.25  for  this  steady 
state,  which  corresponds  to  a  25%  upstream  bias.  The  obvious  location  where  ip  increases  mildly, 
to  Ip  =  0.275,  is  at  the  curvature  discontinuity,  which  correlates  with  the  theoretical  finding  (91). 

Figure  14  b)  presents  the  distribution  of  Mach  nmnber  M,  which  also  reflects  the  exact  solution, 
and  the  corresponding  distributions  of  the  acoustics  and  pressure-gradient  upstream-bias  functions 
a  and  6.  According  to  these  distributions,  the  acoustics  upstream-bias  vanishes  for  M  >  0.7  and 
is  present  with  decreasing  weight  only  in  the  subsonic  region  of  the  flow,  for  M  <  0.7.  For  these 
Mach  numbers,  J  ~  0,  which  corresponds  to  a  centered  approximation  of  the  pressmre  gradient  for 
any  ip.  For  increasing  M  beyond  M  =  0.7,  6  briskly  rises,  which  corresponds  to  a  rapid  upstream- 
bias  growth  in  the  approximation  of  the  pressure  gradient.  For  M  >  1,  S  =  1,  which  makes  the 
pressure-gradient  upstream  bias  weight  equal  to  that  on  the  convection  divergence. 


Figure  14:  Isentropic  Flow:  a)  Energy  and  Controller,  b)  Mach  Number,  a  and  S 
8.2.2  Flow  with  Embedded  Normal  Shock 

A  normal  shock  wave  features  in  this  steady  flow  as  a  result  of  the  subsonic-outlet  pressure  boundary 
condition  p/ptoti„  =  0.84,  imposed  as  an  impulsive  step  decrease  from  the  initial  conditions  and 
for  the  entire  flow  evolution  toward  steady  state.  The  theoretical  solution  places  the  normal  shock 
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at  the  area  ratio  As/j4*in  =  1.09896,  which  corresponds  to  the  interior  of  the  finite  element  with 
node  coordinates  x  =  0.64  and  x  =  0.65,  within  the  computational  domain.  The  exact  shock 
Mach  numbers  are  Mgup  =  1.36989  and  Mgub  =  0.75274,  which  lead  to  the  stagnation  pressure 
and  critical-area  ratios  Ptotout/Ptotin  =  ■^»inM*out  =  0.96537.  The  associated  outlet  area  ratio  and 
Mach  number  are  >loutM»in  =  1-5,  Mout  =  0.45025. 

Figures  15  a)-b)  present  the  convergence  rate  and  steady-state  speed.  The  steady  state  was 
achieved  in  about  25  time  steps,  for  a  total  of  about  50  IRK  cycles,  with  Cmax  =  400  and  a  reduction 
of  the  maximum  residual  by  13  orders  of  magnitude,  down  to  machine  zero.  In  comparison. 
Reference  13  reports  that  a  steady  state  for  the  same  problem  was  achieved  with  a  minimum 
of  175  cycles.  The  calculated  speed  reflects  the  exact  solution,  indicated  with  a  solid  line.  This 
distribution  clearly  shows  the  sonic-point  curvature  discontinuity,  as  well  as  the  normal  shock  which 
is  captured  at  the  correct  location,  within  about  two  elements  with  a  mild  shock-foot  undershoot, 
similar  to  that  in  Reference  13. 


Figure  15:  Shocked  Flow:  a)  Convergence  Rate,  b)  Speed 

The  distributions  of  static  density  and  pressure  in  Figures  16  a)-b)  also  reflect  the  corresponding 
exact  solutions  with  the  curvature  discontinuity  clearly  resolved  and  normal  shock  sharply  captured 
over  only  one  internal  node  with  slight  shock-foot  overshoots.  Significantly,  the  algorithm  has 
correctly  generated  continuous  distributions  for  both  linear  momentum  m  and  enthalpy  H  across 
the  normal  shock,  without  any  shock  bump.  Furthermore,  the  computed  H  remains  accurately 
constant,  as  has  to  be  the  case  for  a  steady  adiabatic  flow. 

Figure  17  a)  presents  the  distributions  of  total  energy  E,  which  visually  coincides  with  the 
exact  solution  in  solid  line,  and  upstream  bias  controller  ip,  with  0.5  <  ^  <  1.0.  The  controller 
remains  essentially  constant  over  smooth  solution  regions,  with  ip  =  0.5,  which  corresponds  to  a 
50%  upstream  bias.  At  the  shock,  ip  rapidly  rises,  reaching  a  ^  ~  0.9  extremum,  which  induces  a 
90%  upstream-bias.  This  increase  in  ip  appears  sudden;  presumably  a  slightly  milder  variation  of  ip 
at  the  shock  could  obviate  the  modest  overshoots  in  p,  p  and  E.  On  the  other  hand,  the  algorithm 
succeeds  in  focusing  an  increased  level  of  upstream-bias,  hence  artificial  dissipation  at  the  shock 
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region  only,  precisely  where  required  for  an  essentially  sharp  and  non-oscillatory  solution. 


Figure  16:  Shocked  Flow:  a)  Density  and  Momentum,  b)  Pressime  and  Enthalpy 


Figure  17:  Shocked  Flow:  a)  Energy  and  Controller,  b)  Mach  Number,  a  and  6 

Figure  17b)  presents  the  distribution  of  Mach  number  M,  which  also  reflects  the  exact  solution, 
and  the  corresponding  distributions  of  the  acoustics  and  pressure-gradient  upstream-bias  functions 
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a  and  5.  According  to  these  distributions,  the  acoustics  upstream-bias  is  present,  with  decreasing 
weight,  only  within  the  inlet  and  outlet  subsonic  regions  of  the  flow.  This  type  of  upstream- 
bias  vanishes  in  the  supersonic  region,  including  the  normal  shock,  hence  it  plays  no  local  role  in 
the  computation  of  the  shock.  Shock  resolution,  therefore,  is  entirely  due  to  the  convection  and 
pressure-gradient  upstream  biases,  which  occur  with  the  same  weight  at  the  supersonic  side  of  the 
shock,  where  5  =  1.  As  M  falls  across  the  shock,  so  does  6,  which  indicates  that  the  upstream 
bias  in  the  approximation  of  the  pressure  gradient  quickly  decreases,  leading  to  an  essentially 
centered  approximation  of  this  gradient  towards  the  outlet.  For  all  the  distributions  computed  for 
this  benchmark,  the  outlet  variations  remain  smooth  and  undistorted;  in  particular  the  calculated 
outlet  pressure  coincides  with  the  imposed  pressure  boundary  conditions,  which  reflects  favorably 
on  the  surface-integral  pressure  enforcement  strategy  delineated  in  Section  6.2. 

8.3  Flows  in  a  Diverging  Nozzle 

These  benchmarks  examine  the  capability  of  the  algorithm  to  calculate  steady  isentropic  and 
shocked  flows  that  involve  a  high-Mach  number  supersonic  region.  The  nozzle  cross-section  area 
distribution  features  a  steep  increase  in  the  diverging  region^^,  as  shown  in  Figure  18,  which  makes 
it  challenging  numerically  to  compute  a  non-oscillatory  shock  located  in  such  a  region. 

The  nozzle  area  ratio  distribution  for  these  benchmarks  is 


=  a  +  b  tanh(  8a:  —  4 ) 

(102) 

with 

0  =  1.39777  , 

b  =  0.34760 

(103) 

and 

=  1.05041  , 

=  1.74514 

A.in 

(104) 

Figure  18:  Variation  of  Area  Ratio  A{x)IA. 
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The  initial  conditions  for  the  gas  correspond  to  an  M  =  1.26  uniform  supersonic  state,  which 
leads  to  the  following  initial  state  throughout  the  nozzle 

/)  =  0.50189  ,  m  =  0.65187  ,  E  =  1.37567  (105) 

The  inlet  flow  is  constrained  supersonic  at  M  =  1.26;  Dirichlet  boundary  conditions  are  thus 
enforced  on  density  p,  linear  momentum  tti  and  total  energy  E.  An  outlet  pressure  boundary 
condition  is  also  imposed  for  the  simulation  of  a  shocked  flow. 

8.3.1  Isentropic  Supersonic  Flow 

For  the  given  initial  conditions,  the  outlet  is  already  supersonic,  hence  no  boundary  conditions  are 
enforced  at  the  outlet.  The  evolution  of  the  flow  from  the  initial  conditions  is  thus  entirely  driven 
by  the  area  source  term  in  the  governing  Euler  equations  (l)-(2). 

Figures  19  a)-b)  present  the  convergence  rate  and  steady-state  speed.  The  acoustics-convection 
upstream  resolution  algorithm  generated  a  steady-state  with  a  reduction  of  the  Tnavimnm  residual 
by  14  orders  of  magnitude,  down  to  machine  zero  in  about  60  time  steps  with  C^ax  =  800.  The 
calculated  speed  reflects  the  exact  solution,  indicated  with  a  solid  line.  This  distribution  clearly 
shows  the  rapid  increase  that  is  triggered  by  the  drastic  nozzle  enlargement. 


Figure  19:  Isentropic  Flow:  a)  Convergence  Rate,  b)  Speed 

The  distributions  of  density,  linear  momentum,  pressure  and  enthalpy  in  Figures  20  a)-b) 
also  mirror  the  corresponding  exact  solutions,  with  swift  expansions  clearly  resolved.  The  algo¬ 
rithm  has  also  correctly  held  constant  the  computed  enthalpy,  which  satisfies  the  steady-adiabatic- 
flow  constant-enthalpy  condition.  Figure  21  a)  presents  the  distributions  of  total  energy  E, 
which  remains  indistinguishable  from  the  exact  solution,  and  upstream  bias  controller  ip,  with 
0.25  <  tp  <  1.0.  Since  the  solution  is  smooth,  %p  stays  virtually  constant  with  ‘ip  ~  0.25  for  this 
steady  state,  which  corresponds  to  a  25%  upstream  bias. 
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Figure  20:  Isentropic  Flow:  a)  Density  and  Momentum,  b)  Pressure  and  Enthalpy 


Figure  21  b)  presents  the  distribution  of  Mach  number  M,  which  also  agrees  with  the  exact 
solution,  and  the  variations  of  acoustics  and  pressure-gradient  upstream-bias  functions  a  and  <5. 
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Figure  21:  Isentropic  Flow:  a)  Energy  and  Controller,  b)  Mach  Number,  a  and  S 

For  a  supersonic  flow,  a  =  0  and  J  =  1.  No  acoustics-matrix  upstream-bias  is  thus  present  in 
this  solution  and  the  entire  flux  divergence  receives  a  uniform  upstream-bias  approximation. 
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8.3.2  Flow  with  Embedded  Normal  Shock 

A  normal  shock  wave  features  in  this  steady  flow  as  a  result  of  the  subsonic-outlet  pressure  boundary 
condition  p/ptoti„  =  0.746,  imposed  as  an  impulsive  step  decrease  from  the  initial  conditions,  for 
the  entire  flow  evolution  toward  steady  state.  The  theoretical  solution  places  the  normal  shock 
at  the  area  ratio  As/A*jn  =  1.35016,  which  corresponds  to  the  interior  of  the  finite  element  with 
node  coordinates  x  =  0.48  and  x  =  0.49,  within  the  computational  domain.  The  exact  shock 
Mach  numbers  are  Mgup  =  1.71319  and  Mgub  =  0.63717,  which  lead  to  the  stagnation  pressure 
and  critical-area  ratios  Ptotout/Ptotin  ==  ■^*in/A»out  =  0.85022.  The  associated  outlet  area  ratio  and 
Mach  number  are  AoutM*in  =  1.74514,  Mout  =  0.43629. 

Figures  22  a)-b)  present  the  convergence  rate  and  steady-state  speed.  The  algorithm  is  certainly 
capable  of  driving  the  maximum  residual  down  to  machine  zero,  with  a  total  residual  reduction 
by  14  orders  of  magnitude,  with  Cmax  =  12;  a  residual  reduction  by  6  orders  of  magnitude  occurs 
within  about  250  time  steps.  Other  initial  conditions,  closer  to  the  steady  state,  would  presumably 
lead  to  faster  convergence.  These  results,  however,  compare  favorably  with  those  in  Reference  13, 
which  reports  lack  of  convergence  with  Roe’s  algorithm  for  this  benchmark  and  with  a  similar 
computational  solution.  This  comparison  rests  on  the  observation  that  the  acoustics-convection 
upstream  resolution  algorithm  is  not  a  purely  flux- vector  splitting  algorithm,  but  it  also  uses,  like 
Roe’s  algorithm,  a  Riemann  solver,  ‘though  applied  to  the  acoustics  equations.  The  calculated 
speed  reflects  the  exact  solution,  indicated  with  a  solid  line,  and  clearly  shows  the  expected  rapid 
rise  preceding  the  shock  as  well  as  an  excellent  calculated  normal  shock,  captured  over  only  one 
node. 


Figure  22:  Shocked  Flow:  a)  Convergence  Rate,  b)  Speed 

The  distributions  of  static  density  and  pressure  in  Figures  23  a)-b)  also  reflect  the  corresponding 
exact  solutions,  with  rapid  changes  in  these  two  variables  clearly  resolved  and  normal  shock  again 
sharply  captured  over  one  internal  node.  In  harmony  with  the  previous  benchmark  results,  also 
for  this  problem  has  the  algorithm  correctly  generated  continuous  distributions  for  both  linear 
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momentum  m  and  enthalpy  H  across  the  normal  shock,  without  any  shock  bump.  Furthermore, 
the  computed  H  remains  again  constant,  as  has  to  be  the  case  for  a  steady  adiabatic  flow. 

Figure  24  a)  presents  the  distributions  of  total  energy  E,  which  visually  coincides  with  the 
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Figure  23:  Shocked  Flow:  a)  Density  and  Momentum,  b)  Pressure  and  Enthalpy 


Figure  24:  Shocked  Flow:  a)  Energy  and  Controller,  b)  Mach  Number,  a  and  6 
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exact  solution  in  solid  line,  and  upstream  bias  controller  if),  with  0.4  <  i/i  <  1.5.  The  controller 
remains  essentially  constant  over  smooth  solution  regions,  with  V’  =  0.4,  which  corresponds  to 
a  40%  upstream  bias.  At  the  shock,  ip  rapidly  rises  and  reaches  a,  ip  0.96  extremum,  which 
induces  a  96%  upstream-bias.  This  increase  in  ip  appears  less  sudden  at  its  inception  than  that  in 
Figure  17  a),  which  reinforces  the  conjecture  that  a  relatively  smooth  increase  of  V*  at  a  normal 
shock  can  obviate  modest  shock  overshoots.  No  overshoots  are  present  in  this  solution,  and  as 
Figure  24  a)  bears  out  the  algorithm  succeeds  in  focusing  an  increased  level  of  upstream-bias, 
hence  artificial  dissipation,  at  the  shock  region  only,  precisely  where  required  for  an  essentially 
sharp  and  non-oscillatory  solution.  Figure  24  b)  presents  the  distribution  of  Mach  number  M, 
which  also  reflects  the  exact  solution,  and  the  corresponding  distributions  of  the  acoustics  and 
pressure-gradient  upstream-bias  functions  a  and  6.  According  to  these  distributions,  the  acoustics 
upstream-bias  is  present  only  in  the  outlet  subsonic  regions  of  the  flow.  This  type  of  upstream- 
bias  vanishes  in  the  supersonic  region,  including  the  normal  shock,  hence  it  plays  no  local  role 
in  the  computation  of  the  shock.  Also  in  this  case,  therefore,  is  shock  resolution  entirely  due 
to  the  convection  and  pressure-gradient  upstream  biases,  which  occur  with  the  same  weight  at 
the  supersonic  side  of  the  shock,  where  5  =  1.  As  M  falls  across  the  shock,  so  does  5,  which 
indicates  that  the  upstream  bias  in  the  approximation  of  the  pressure  gradient  quickly  decreases, 
leading  to  an  essentially  centered  approximation  of  this  gradient  towards  the  outlet.  Also  for  all 
the  distributions  computed  for  this  benchmark  are  the  outlet  variations  smooth  and  undistorted;  in 
particular  the  calculated  outlet  pressure  coincides  with  the  imposed  pressure  boundary  conditions, 
which  again  reflects  favorably  on  the  surface-integral  pressure  enforcement  strategy. 

9  Concluding  Remarks 

The  acoustics-convection  upstream  resolution  algorithm  rests  on  the  physics  and  mathematics  of 
acoustics  and  convection.  It  introduces  a  decomposition  of  the  flux  vector  Jacobian  into  acoustics 
and  convection  matrix  components,  for  general  equilibrium  equations  of  state,  and  generates  the 
upstream  bias  at  the  differential  equation  level  before  any  discrete  approximation.  This  formulation 
generates  a  characteristics-bias  flux  which  generalizes  in  the  continuum  the  traditional  upwind- 
scheme  numerical  fluxes. 

A  natural  finite  element  discretization  of  the  characteristics-bias  flux  generates  an  essentially 
centered  approximation  of  the  Euler  flux  divergence,  in  the  form  of  a  non-linear  combination  of  up¬ 
stream  diffusive  and  downstream  anti-diffusive  flux  differences,  with  greater  bias  on  the  upstream 
diffusive  flux  difference.  The  algorithm  induces  this  upstream  bias  and  associated  artificial  diffu¬ 
sion  mostly  locally  in  regions  of  solution  discontinuities,  whereas  it  decreases  the  upstream-bias  in 
regions  of  solution  smoothness.  This  solution-driven  variable  upstream-bias  thus  minimizes  artifi¬ 
cial  dissipation  via  an  upstream-bias  controller  that  depends  on  the  jumps  of  solution  slopes.  The 
study  in  this  paper  has  implemented  the  algorithm  using  a  linear  approximation  of  fluxes  within 
two-noded  cells,  without  any  MUSCL-type  local  extrapolation  of  variables. 

The  computational  results  have  validated  accuracy  and  essential-monotonicity  performance 
of  the  acoustics-convection  upstream  resolution  algorithm  for  transient  and  steady  smooth  and 
shocked  flows.  The  algorithm  generates  essentially  non-oscillatory  solutions  and  automatically 
preserves  a  constant  enthalpy  as  well  as  smoothness  of  both  enthalpy  and  mass  flux  across  steady 
normal  shocks.  These  results  have  also  validated  an  intrinsically  stable  pressure  boundary  condition 
procedure  at  a  subsonic  outlet.  This  procedure  directly  enforces  the  outlet  pressure  within  the 
surface  integral  that  emerges  in  the  momentum-equation  weak  statement.  The  computed  solutions 
at  nozzle  outlets  remain  smooth  and  undistorted  and  mirror  the  exact  reference  solutions,  which 
bears  out  the  reliability  of  this  pressure  boimdary  condition. 
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According  to  the  solution-driven  numerical  values  of  the  upstream-bias  controller,  the  computed 
smooth  and  shocked  solutions  resulted  from  a  mostly  centered  discretization.  This  finding  indicates 
that  a  uniformly  fully  upwind  formulation  is  not  strictly  necessary  within  a  characteristics-bias 
algorithm  to  generate  essentially  sharp  and  non-oscillatory  solutions.  The  algorithm,  therefore, 
succeeds  in  both  reducing  artificial  dissipation  in  regions  of  smooth  flow,  for  higher  accuracy,  and 
focusing  an  increased  level  of  upstream-bias,  hence  artificial  dissipation,  at  the  shock  regions  only, 
precisely  where  required  for  stability  and  sharp  shock  capturing. 

The  computational  results  agree  with  the  reference  solutions,  with  curvature  discontinuities 
exactly  calculated  and  sharp  normal  shocks  captured  over  one  or  two  points.  The  finite  ele¬ 
ment  acoustics-convection  upstream  resolution  algorithm,  therefore,  delivers  solutions  that  are  as 
sharp  and  non-oscillatory  as  those  generated  by  local-extrapolation  flux-vector  and  flux-difference 
splitting  schemes^^.  This  characteristics-bias  algorithm,  however,  admits  a  straightforward  im¬ 
plicit  implementation,  features  a  computational  simplicity  that  parallels  a  traditional  centered 
discretization,  and  rationally  decreases  superfluous  artificial  diffusion.  Ongoing  work  is  complet¬ 
ing  an  intrinsically  multi-dimensional  and  infinite-directional  acoustics-convection  algorithm,  to  be 
detailed  in  a  future  paper. 
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